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drifters'.  To  analyze  the  irregular  spaced  drifter  tracks,  objective  analysis  techniques  were  applied  to  produce 
optimal  estimations  of  surface  current  fields  on  a  regularly  spaced  grid.  The -transformation  of  a  data  set  onto 
a  regularly  spaced  grid  allows  for  the  use  of  many  other  standard  computer  analysis  programs.  A  computer 
program  using  this  technique  has  previously  been  successfully  applied  to  large  oceanographic  data  sets  of  many 
drifter  tracks.  The  objective  analysis  program  was  converted  for  use  on  Hewlett  Packard  microcomputers 
and  applied  to  a  much  more  limited  data  set  from. six  drifters.  The  results  being  that  objective  analysis  can 
effectively  work  with  small  data  sets. 

The  effect  of  very  limited  data  sets  on  the  outcome  of  the  estimated  fields  was  checked  by  removing  drifters, 
one  at  a  time,  from  the  data  set  of  sii/ rlrifters.  The  fields  generated  by  the  reduced  number  of  drifters  were 
then  compared  to  the  results  with  the  'best  estimate"  field  determined  by  all  six.  drifters.  The  results  from 
this  exercise  are  that  the"  placement  of  the  buoys  or  drifters  relative  to  the  flow  field  features  Will' greatly 
influence  how  well  the  total  flow  field  is  estimated.  A  few  well  placed  buoys  wiH  provide  better  information 
than  many  poorly  placed  buoys.  Therefore,  the  use  of  remotely  sensed  data  will  aid  in  the  determination  of  the 
optimal  placement  of  buoys.  •  , .-  i  ,  , 
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Chapter  1 
INTRODUCTION 


1.1  BACKGROUND 

The  Oceanography  Branch  has,  as  part  of  the  Improvement  of 
Probability  of  Detection  in  Search  and  Rescue  Project,  an 
ongoing  investigation  into  the  drift  and  leeway  of  survivors  and 
survivor  crafts  (see  Paskausky,  1986).  Murphy  et  ad  1982; 
Anderson,  1984;  and  Murphy  and  Allen,  1985;  evaluated  the  Coast 
Guard's  operational  search  planning  computer  models.  They 
showed  that  model  predictions  of  survivors  and  survivors  craft's 
drift,  based  on  either  historical  current  files  or  surface 
current  generated  by  large  scale  winds  from  U.S.  Navy  Fleet 
Numerical  Oceanography  Center  (FNOC) ,  rarely  predicted  the 
actual  drift  of  survivor  craft  as  simulated  by  freely  drifting 
buoys.  At  present,  a  single  radio-direction  finder  type  datum 
marker  buoy  (RDF-DMB)  is  used  to  measure  surface  currents 
during  a  search.  The  initial  position  of  the  buoy  is  established 
when  the  RDF-DMB  is  dropped  from  an  aircraft.  However,  to 
obtain  a  velocity  datum,  a  second  position,  at  a  later  time, 
must  be  determined  for  the  RDF-DMB.  This  requires  the  aircraft 
to  break  off  the  search  for  the  survivors  and  search  for  the 
RDF-DMB,  causing  the  loss  of  valuable  time  and  fuel.  The 
important  velocity  measurements  necessary  for  an  accurate 
prediction  of  the  drift  of  the  survivors  are  therefore  limited 
to  one  or  two  points.  The  development  of  VHF  Loran-C  buoys 
(Allen,  Eynon,  and  Robe,  1987) ,  demonstrate  that  coded  buoys 
reporting  their  Loran-C  positions  every  30  minutes  will  deliver 
far  more  positions,  more  accurately  and  therefore  more  surface 
current  data  to  the  search  planner.  This  provides  the  search 
planner  a  clearer  picture  of  the  surface  currents.  The 
increased  data  rates  requires  a  more  involved  analysis  than  the 
data  from  the  RDF-DMB' s.  This  report  takes  a  first  look  at  an 
analysis  technique  for  freely  drifting  buoys  and  other  sources 
of  randomly-spaced,  surface  current  data.  This  technique  is 
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known  as  objective  analysis. 


The  use  of  statistical  models  for  the  optimal  estimation  of 
oceanic  fields  '•.ven  limited  measurements  is  known  as  objective 
analysis  (Bretherton,  et  al,  1976) .  This  technique  has  been 
applied  by  Gandin  (1965)  to  analyze  the  wind  and  pressure  field 
in  the  atmosphere,  and  used  routinely  for  the  preparation  of 
numerical  weather  prediction.  It  is  a  very  important  tool  in 
oceanography  for  both  analysis  and  for  observational  array 
design  (White  and  Bernstein,  1979?  Clancy,  1983  and  Robinson  and 
Leslie,  1985) . 

The  models  described  here  are  the  extension  of  space-time 
objective  analyses,  and  serve  as  the  statistical  component  of 
the  Harvard  Ocean  Descriptive  Predictive  System  (Robinson  and 
Leslie,  1985) .  They  can  be  used  as  an  interpolation  scheme  both 
spatially  and  temporally  to  provide  initial  and/or  boundary 
conditions  for  the  dynamic  model  (Tu,  1981) ,  or  used  to  forecast 
the  oceanic  fields  (Carter,  1983) .  The  essential  assumption  in 
practical  use  is  sufficient  statistical  knowledge  about  how  the 
fields  are  related  both  in  time  and  space. 

In  this  study,  objective  analysis  was  utilized  to 
nowcast/ forecast  the  surface  currents  using  limited  drift 
measurements.  Several  types  of  freely  drifting  buoys  that 
transmit  their  positions  to  a  shore-based  receiver  have  been 
used  by  the  Oceanography  Branch.  A  local  positioning  system 
using  the  Microwave  Tracking  System  (MTS)  was  used  with  surface 
drifters  in  leeway  studies.  Both  polar-orbiting,  NOAA/TIROS 
series  satellite  and  Loran-C  positioning  were  used  with  surface 
drogued  buoys  during  drift  studies  on  the  continental  shelf. 
All  of  these  studies  (Murphy  and  Allen,  1985;  Allen,  Eynon, 
and  Robe,  1987;  Nash  and  Willcox,  1988)  typically  contained  data 
sets  in  the  form  of  a  time  series  of  buoy  positions  for  one  to 
six  drifters.  A  time  series  of  a  freely-drifting  buoy's 
positions  is  called  a  drift  track.  The  purpose  of  this  study 
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was  to  apply  the  technique  of  objective  analysis  to  one  of  our 
data  sets.  A  computer  program  using  this  technique  has  been 
successfully  applied  to  large  oceanographic  data  sets  of  drifter 
tracks  (Robinson  and  Leslie,  1985) .  These  large  data  sets, 
however,  included  many  more  drifter  tracks  than  are  typically 
used  in  our  studies.  Therefore,  the  objective  analysis  program 
was  converted  for  use  on  Hewlett  Packard  microcomputers  and 
applied  to  a  data  set  from  an  MTS  leeway  study. 

1.2  BASIC  THEORY 

The  basic  equations  of  the  objective  analysis  procedure  are  well 
described  in  the  paper  by  Bretherton,  et  al.  (1976),  to  which 
readers  are  referred  for  the  details.  In  this  section  the 
special  case  of  scalar  objective  analysis  is  described  and  the 
general  matrix  forms,  which  is  applicable  to  analysis  of  both 
scalar  and  vector  fields,  are  presented. 

In  scalar  objective  analysis,  the  value  ©p  of  a  scalar  variable 
6(x,y)  at  a  general  point  Xp  (x,y)  is  estimated  from  measurements 
at  data  points  (x,y) ,  where  m  *  1,  ...,  N.  A  basic 

assumption  of  the  procedure  is  that  the  field  6  is  stationary  and 
homogeneous  and  has  known  mean  and  covariance.  The  field 
statistics  are  typically  estimated  from  a  combination  of  a  priori 
assumptions  and  a  sample  of  measurements. 

Each  measured  value  at  point  Xjjj  is  considered  to  be  the  true 
value  0^  of  the  field  at  that  point  plus  some  random  noise  em, 
due  to  sampling  or  instrumental  error: 

V9m+,V  . . Nl  '  (1) 

where  N  is  the  number  of  measurements.  The  errors  e_  are  assumed 

m 

to  be  uncorrelat  ’  with  the  field  8  and  with  one  another  and  have 
known  variance  ? 
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em  n 


0 


(2) 


em  e  n  £  8  mu  /(m, n=l,...,N),  (3) 

where  the  f unctionS  „  has  value  of  unit  when  m  equals  n,  but  zero 

mn 

otherwise.  What  the  above  assumes  is  that  there  are  no 
systematic  or  calibration  errors. 


To  reconstruct  the  entire  field  the  estimator  ep  for  0(x, y)  at 

each  point  is  found.  The  Gauss-Markov  theorem  states  that  the 
P 

optimal  linear  estimator  for  0(x,y)  at  each  position  Xp(x,y)  is: 


where  a“L  are  elements  of  the  inverse  of  the  matrix  of 
mn 

covariances  between  all  NxN  pairs  of  observations,  and  C  is  the 

pm 

covariance  between  the  field  value  ©_  to  be  estimated  and  the  nth 

P 

measurement.  The  covariances  are  generally  assumed  to  be  a 
function  of  the  distance  between  the  field  positions. 
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c  -0  4>  =  F  (x  -  x  ) 

pm  pm  pm 


(6) 


For  given  estimation  and  observation  positions,  A_  and  c  are 

mn  pm 

constants.  Equation  (4)  therefore  expressed  the  estimate  for  the 
field  variable  as  a  linear  combination  of  the  measurements, 
weighted  by  the  covariances. 

The  minimum  mean  square  error  for  the  estimator  is  the  variance 

of  the  error  in  §  given  by: 
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The  first  term  is  the  overall  natural  variation  in  the  field  ©, 
while  the  second  term  is  a  measure  of  information  provided  by  the 
measurements . 

A  more  general  vector  form  for  equations  (4)  and  (7)  is: 


A 

6 


-l 

c  $0  c  $  $ 


(8) 


Ce  -  C0  -  , 


(9) 


where  £  is  the  estimated  vector  describing  the  entire  field.  C^Q 
is  the  cross-covariance  matrix  between  the  field  vector  and  the 
observation  vector.  is  the  auto-covariance  matrix  for  the 
vector  of  observations  0.  CQ  is  the  natural  covariation  of  the 
field  as  defined  by  the  a  priori  covariance  function.  is  the 
error  variance.  The  superscript  T  denotes  the  matrix  transpose. 

Equations  (8)  and  (9)  are  general  and  can  be  applied  to  a  two- 
dimensional  or  higher  vector  variable.  For  example,  given  two 
component  (u,v)  observation  of  velocity  at  N  points,  the 
covariance  matrix  now  includes  the  auto-covariances  of  each 
observational  component  and  cross-covariances  between  the  two 
components.  The  matrix  C^Q  includes  the  covariance  between  each 
observation  component  and  each  component  of  the  estimated  field 
variable.  Vectors  £  and  0  include  values  for  each  component  at 
each  estimation  and  observation  point,  respectively.  It  should 
be  pointed  out  that  the  theorem  does  allow  replacing  the 
covariance  by  the  correlation  functions  in  practical  use. 
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2 . 1  INTRODUCTION 


Chapter  2 
METHODS 


The  covarinace  matrices  required  in  the  solution  of  the 
objective  analysis  as  described  in  the  Basic  Theory  section  are 
defined  by  correlation  functions.  The  correlation  function  is 
the  planar  (x,y)  and  temporal  (t)  distribution  of  the  weighting 
factors  for  estimating  an  unknown  value  from  the  surrounding 
observations.  Essentially,  closer  observations  are  weighted 
more  in  determining  an  unknown  value  than  more  distance  values. 
At  sufficiently  large  distances  and  time  lags  the  observations 
will  have  no  appreciable  influence  on  the  interpolated  value. 
Thus,  a  limit  over  which  the  observations  are  not  included  in 
the  interpolation  is  assumed.  This  greatly  simplifies  the 
computation.  The  correlation  function  can  be  determined 
directly  from  a  data  set  that  is  sufficiently  large.  However, 
when  the  data  set  is  small,  the  mathematical  form  of  the 
correlation  function  is  artificially  imposed. 

2.2  THE  DETERMINATION  OF  CORRELATION  FUNCTION 

The  essential  assumption  in  objective  analysis  is  that  the 
correlation  functions  are  known.  These  functions  can  be 
calculated  from  the  measurements  by  statistical  approaches,  if 
there  are  large  and  intensive  data  sets  available.  The 
determination  of  these  correlation  functions  are  very  difficult 
when  the  observed  data  are  sparse.  An  alternate  approach  is  to 
use  a  fitted  analytical  formula  which  represents  the  correlation 
function  of  the  limited  data.  The  correlation  function  should 
satisfy  the  following  two  requirements:  (1)  symmetric  to  the 
lag,  r,  i.e.,  C(r)=C(-r),  where  r  =  (delta_x,delta_y,delta_t)  is 
the  space-time  lag,  and  (2)  in  positive  definite  form.  It  is 
important  to  note  that  the  calculated  correlation  function  may 
contain  some  small  scale  noise  which  must  be  filtered.  This 
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smoothed  correlation  function  is  then  used  as  a  "look-up  table" 
for  the  objective  analysis. 

2.3  THE  ELIMINATION  OF  DISTANT  DATA 

The  correlation  function  drops  off  with  increasing  space  and 
time  lags.  In  other  words,  the  distant  (in  space  and  time) 
observations  have  very  little  influence  on  an  interpolation 
point  when  compared  to  nearby  points.  So  the  elimination  of  the 
distant  observations  could  make  the  computation  very  efficient. 
Usually,  limiting  the  domain  of  influence  (in  space  and  time)  of 
an  interpolation  point  is  sufficient  to  reduce  the  matrices  used 
in  the  analysis  to  a  reasonable  size,  but  with  a  large  data  set 
the  limited  influence  domain  may  still  contain  a  lot  of  data. 
For  the  large  data  case,  the  maximum  influence  points  are 
"limited"  to  less  than  th*5  full  number  of  data  points  in  the 
domain.  The  data  points  which  have  the  highest  correlation  with 
the  interpolation  points  or  the  points  that  individually  give 
the  lowest  error  estimate  are  ordered  by  rank,  and  then  the  top 
"limited"  values  are  used  in  the  calculation.  These  are  the 
optimal  data  points. 

2.4  THE  CORRELATION  FUNCTION  USED  IN  THIS  STUDY 

The  correlation  function  used  in  this  present  study  is: 

-R 

C  (R)  =  e  2  2  ? 

where  R=  V  f(_Ax  )  +  (  A^_)  +  (_Ajt  )2  1  (10) 

LVxfold/  vyfold/  vtfold/  J 

where  delta_x,  delta_y,  and  delta_t  are  the  components  of  the 
difference  from  the  data  point  to  the  interpolated  point;  and 
where  xfold,  yfold,  and  tfold  are  the  e-folding  scales  of  the 
correlation  function. 

Equation  (10)  for  x  and  y  is  shown  in  Figure  1,  for  values  of 
xfold,  and  yfold  equal  1.  At  x,y  =  1,0  and  x,y  =  0,1  the 
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correlation  (z)  value  is  0.368,  this  is  called  the  e-folding 
scale.  The  e-folding  scale  is  shown  in  Figure  1  as  the  circle 
at  z  =  0.368.  As  the  correlation  function  approaches  zero 
asymptotically,  the  outer  portion  is  cut  off  by  the  limit  set  in 
Equation  (6)  .  The  appropriate  values  for  the  e-folding  scales 
are  the  physical  scales  of  the  motions.  This  is  the  scale  in 
time  and  space  of  the  coherence  of  the  dominant  motions.  The 
e-folding  scale  is  symmetrical  about  its  axis,  but  does  not 
have  to  be  equal  about  all  axes.  Figure  2  is  the  correlation 
function  where  xfold  =  1  and  yfold  =  3 .  In  Figure  3  the  xfold  = 
3  and  yfold  =1,  in  these  cases  the  e-folding  scale  are  ellipses 
at  z  =  0.368.  These  simple  examples  can  be  transferred  to  a 
more  realistic  case,  e.g.,  the  field  in  question  is  alongshore 
velocity  on  the  continental  shelf,  the  coherence  is  very  long 
(-200km)  in  the  along  shore  (y)  direction  and  short  in  the  cross 
shelf  (x)  direction  (~20kms)  and  about  two  days  in  (t)  time. 
These  values  could  be  used  as  a  first  guess  of  the  e-folding 
scales  for  x,  y  and  t.  The  appropriate  values  are  based  on  the 
physical  scales  of  the  value  used. 

2 . 5  PROGRAM  DESCRIPTIONS 

The  SOA/VOA  (Scalar/Vector  Objective  Analysis)  Programs  computer 
software  package  consists  of  twelve  programs  and  related 
documentation,  which  were  originally  written  by  Carter  (1983)  in 
FORTRAN  language  for  the  VAX  series  mainframe  computers.  They 
were  rewritten  by  Dr.  Sun,  in  Hewlett  Packard  (HP)  BASIC  and 
designed  to  run  on  Hewlett  Packard  9000  series  200  mini¬ 
computers.  Many  enhancements  were  added  to  make  the  programs 
easier  to  use  (for  example,  better  prompts  for  inputting  the 
parameters  and  outputting  the  results,  easier  error  checking  on 
inputs,  etc.).  The  distribution  disc  labeled  "DISC_1"  contains 
the  twelve  programs  and  related  documentation. 

The  package  contains  programs  for  displaying  the  original  data, 
setting  up  the  objective  analysis,  both  a  scalar  and  a  vector 
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analysis  program,  and  programs  for  displaying  and  analyzing  the 
results  of  the  objective  analysis.  There  are  three  programs 
that  display  the  data  from  the  drifters.  The  u  and  v  components 
of  velocity  can  be  plotted  as  a  time  series,  or  individually  on 
the  x,y  plane.  The  PREPARE  and  GET_DATA  programs  set  up  the 
drifter  data  for  input  into  the  main  objective  analysis 
programs.  The  GET_IP  programs  converts  the  original  coordinate 
system  to  the  interpolated  positions  of  the  array  used  in  the 
objective  analysis.  Either  of  two  programs  provide  the 
correlation  function  for  the  objective  analysis.  COM_COR 
calculates  the  correlation  function  from  the  observed  data  set 
and  the  subroutine  FNCOV_F  calculates  the  function  from  a  given 
mathematical  formula.  Either  the  Scalar  Objective  Analysis 
(SOA)  or  the  Vector  Objective  Analysis  (VOA)  are  the  main 
programs  for  determining  the  interpolated  fields.  Working 
programs  are  stored  under  the  names  SOA_W  and  VOA_W.  Two 
programs  are  for  displaying  the  results  of  the  objective 
analysis.  PRINT_OA  shows  the  individual  numeric  values  on  the 
CRT  and  prints  them  on  the  printer.  PLOT_OA  draws  the  contours 
of  the  results  and  displays  them  on  the  CRT  or  prints  then  on  a 
printer  or  plotter.  The  difference  between  two  fields  generated 
by  separate  runs  of  the  objective  analysis  is  calculated  by  the 
program  SOA_DIFF. 

A  brief  description  of  the  programs  follow. 

The  three  programs  that  display  the  data  from  the  drifters  are: 
LINDAT_UV  PROGRAM: 

plots  the  u-  and  v-component  of  the  current  speeds  as  a  function 
of  time.  The  u-  and  v-speed  are  indicated  as  solid  and  dotted 
lines  respectively. 

PLOT_PUS  PROGRAM: 

plots  the  drifter  position  and  the  u-component  of  the  current 
speed  at  a  specific  time. 
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PLOT  PVS  PROGRAM: 


plots  the  drifter  position  and  the  v-component  of  the  current 
speed  at  a  specific  time. 


The  three  programs  that  set  up  the  drifter  data  and  the 
interpolated  array  for  input  into  the  main  objective  analysis 
programs  are: 

PREPARE  PROGRAM: 


prepares  the  observation  data  for  test  run  by  converting  the 
original  data  (x,y,t)  to  data  files  by  drifter  of  x,y,t,u,  and 
v. 

GET_DATA  PROGRAM: 

reads  in  the  observed  data,  and  sets  up  the  input  and  output 
files  for  the  objective  analysis  program. 

GET_IP  PROGRAM: 

converts  the  original  coordinate  system  (e.g.  Florida  state 
plane)  to  inter/extrapolation  positions  of  the  objective 
analysis  array. 

Two  programs  provide  the  correlation  function  for  the  objective 
analysis,  COM_COR  program  and  FNCov_t  subroutine  or  by  the 
FNCov_f  subroutine. 

COM_COR  PROGRAM: 

calculates  the  correlation  function  from  the  observed  data,  and 
stores  them  in  the  form  of  a  "look-up"  table. 
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FNCov  t  SUBROUTINE: 


gets  the  correlation  function  from  the  "look-up"  table. 


FNCov_f  SUBROUTINE: 

calculates  the  correlation  function  from  a  fitted  analytical 
formula. 

SOA/VOA  PROGRAMS : 

The  "SOA"  and  "VOA"  programs  play  the  key  roles  in  the  scalar 
and  vector  objective  analyses  respectively.  Each  program  is 
composed  of  one  main  program  and  several  sub-programs.  The  main 
program  prompts  for  inputting  parameters  from  keyboard  and  reads 
in  the  measurements  and  inter/extrapolation  position  data,  calls 
sub-programs  to  do  the  objective  analysis,  then  finally  outputs 
the  estimated  fields  and  expected  errors  to  disk  in  Drive  No.  1. 

SOA  W  PROGRAM: 


The  working  program  of  the  scalar  objective  analysis. 

VOA_W  PROGRAM: 

The  working  program  of  the  vector  objective  analysis. 

Two  programs  display  either  the  individual  numerical  results  or 
the  contours  of  results. 

PRINT_OA  PROGRAM: 

displays  the  individual  results  of  the  objective  analysis 
results  on  the  screen.  It  also  outputs  the  results  to  the 
printer,  if  desired. 
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PLOT  OA  PROGRAM: 


contours  the  objective  analysis  results  and  then  plots  to  the 
screen.  It  also  dumps  the  graphics  to  the  printer  or  plotter, 
if  desired. 

SOA  DIFF  PROGRAM: 


calculates  the  difference  between  two  15  by  15  arrays  (File  1  - 
File  2)  and  determines  the  Cumulative  Measured  Error  Index 
(Equation  9) . 

A  brief  description  of  the  major  subprograms  of  the  SOA/VOA 
program  follow. 

DIAG  SUBPROGRAM: 

calculates  the  statistical  parameters  of  inputs  such  as  the 
mean,  variance,  root  mean  square,  minimum  and  maximum. 

EST_MEAN  SUBPROGRAM: 

calculates  the  estimated  mean,  and  then  removes  the  estimated 
mean  from  the  input  array. 

FNIER  SUBPROGRAM: 

returns  error  code  for  subprogram  INVMTX. 

GET  RD  SUBPROGRAM: 


gets  the  data  before  and  after  a  given  time  and  also  within  a 
given  spatial  radius  from  the  domain  reference  point. 

INVMTX  SUBPROGRAM: 
inverts  the  Matrix. 
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SCALAR  OA  SUBPROGRAM: 


The  scalar  space-time  objective  analysis  routine. 

SELECT  SUBPROGRAM: 

eliminates  the  distant  (in  space  and  time)  data  points  and 
selects  the  nearest  data  points  to  an  inter/extrapolation  point 
X, Y,  and  T.  The  number  of  optimal  data  points  is  restricted  to 
the  maximum  number  of  influential  points. 

SET_INVA  SUBPROGRAM: 

sets  up  the  correlation  function  for  the  input  data  given  the 
optimal  positions  and  times,  it  returns  the  inverted  correlation 
function  matrix. 

SORT  SUBPROGRAM: 

sorts  the  index  and  correlation  in  descending  order. 

VECTOR  OA  SUBPROGRAM: 


The  vector  space-time  objective  analysis  routine. 

2.6  USER'S  INSTRUCTIONS  AND  EXAMPLES 

In  this  section  the  reader  is  shown  how  to  min  the  programs  by 
using  the  realistic  field  data. 

The  procedures  to  run  the  SOA  programs  are  as  follows:  (note: 
hit  <ENTER>  key,  after  you  answer  the  program  prompts.) 

Step  1 :  Boot  the  system. 

Step  2:  Insert  the  distribution  disk  labeled  DISC_1  into  Drive 

No.  0. 
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Step  3:  Insert  the  data  disk  labeled  DISC_2  into  Drive  No.l. 

Step  4;  Prepare  the  observational  data. 

The  drifter  position  data  named  "LDIOIAPR"  in  DISC_2  is  used 
here  as  an  example.  The  "LDIOIAPR"  contains  positions  of  ten 
drifters  and  four  boats  as  a  function  of  time.  The  program 
PREPARE  calculates  the  drifting  u-  and  v-speed  for  each  drifter 
and  boat,  and  stores  the  drifter  positions,  times  and  calculated 
current  components  by  using  the  format  of  "X,Y,T,U,V".  Files 
Drifter_l  to  Drifter_10  are  the  drifter  data,  while  Files 
Drifter_ll,  12,  13,  and  14  are  the  drifting  speed  data  for  boats 
Little  Lady,  Whaler,  Joanna,  and  R/V  Oceaneer  IV  (drifter) 
respectively . 

The  inter/extrapolation  of  the  u-speed  from  the  drifter  data  is 
the  objective.  The  program  GET_IP  calculates  the  area  of 
interest  which  is  a  square  domain  with  15x15  grid  points.  The 
inter/ extrapolation  position  data  are  saved  in  file  IP_POS.  The 
correlation  function  used  in  our  present  study  is  equation  (10) . 

The  first  test  run  used  data  from  6  drifters  including  Drifters 
No.  5,  6,  7,  8,  9,  and  10.  The  GET_DATA  program  was  used  to  get 
these  drifter  data  and  store  them  in  file  "OBSDATA6" .  Different 
runs  can  be  made  by  changing  the  file  names  in  the  "CREATE  BDAT" 
lines  1042,  1043,  1044,  and  the  "ASSIGN  @Pathout"  line  1047  and 

commenting  out  the  appropriate  drifter  "DATA"  lines  (1015-1040) 
and  changing  the  upper  limit  of  "N"  in  line  1050  (FOR  N  =  1  to 
n)  . 

Step  5:  Get  a  working  program. 

1.  Insert  the  distribution  disk  labeled  DISC_1  into  Drive  no.  l 

2.  Load  "SOA"  program.  Type  LOAD  "SOA”,  hit  <ENTER>  key. 
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3.  Load  "FNCov_f"  or  "FNCov_t"  subprograms.  Type  LOADSUB 

"FNCov_f " (or  FNCov_t") ,  hit  <ENTER>  key. 

4.  Make  correction  on  line  4655  (ENTER  @Path_in;  X,  Y,  T,  Phi, 

Dummy)  for  observed  data  format,  and  line  4775  (ENTER 
@Path_in;  X,  Y)  for  inter/extrapolation  position  data 
format,  if  needed. 

Step  6;  Store  this  working  program  for  future  use,  if  desired. 
Type  STORE  "SOA_W",  hit  <ENTER>  key. 

Step  7 :  Now,  you  are  in  a  position  to  run  the  Scalar  Objective 

Analysis  and  interested  in  inter/extrapolation  the 
field  on  14:30:00,  1  April  1S85.  Hit  <ENTER>/<RUN> 

key  to  start,  the  program  will  prompt 

"Do  you  need  the  documentation?  (answer  Y/N,  for  yes/no) " 

Type  "N"  for  no  documentation,  and  begin  the  computation. 

1.  "Enter  time  the  analysis  to  make.  (DD  MMM  YY,HH:MM:SS) " 

Type  1  APR  1985,14:30:00 

2.  "Enter  number  of  objective  analyses  to  make." 

Type  1 

3.  "Enter  time  interval  for  extrapolation  in  time." 

Type  0 

4.  "Enter  the  maximum  distance  radius  from  the  reference  point 
of  domain." 

Type  1E+6 

5.  "Enter  the  maximum  time  radius  before  and  after  the  time  of 
the  analysis." 

Type  0 
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6.  "Enter  the  maximum  space  and  time  lags." 

Type  5E+4 , 0 

7.  "Enter  the  maximum  number  of  influential  points." 

Type  6 

8.  "Enter  the  X  direction  e-folding  Scale." 

Type  1E5 

9.  "Enter  the  Y  direction  e-folding  Scale." 

Type  1E5 

10.  "Enter  the  Time  e-folding  Scale." 

Type  1000 

11.  "Enter  the  observed  data  file  including  file  specifier." 

Type  OBSDATA6 :  ,  7  0  0 , 1 

12.  "Enter  the  interpolation  position  file  including  file 
specifier. " 

Type  IP_POS: ,700, 1 

13.  "Enter  the  file  specifier  for  OA  forecast  fields." 

Type  S0A_FCST6U : 700 , 1 

14.  "Enter  the  file  specifier  for  OA  error  variance  fields." 

Type  SOA_EVAR6U : , 7  0  0 , 1 

Step  8:  The  working  program  "SOA_W"  will  output  the  forecast 
field  and  expected  errors  fields  to  Disc_2  in  Drive  #1. 

Step  9:  Output  the  results.  Use  programs  PLOT_OA  and  PRINT_OA 
to  get  the  hard  copy  of  the  analysis  results. 

Repeat  steps  5-9  to  calculate  the  v-component  speed,  after 
making  changes  to  line  4655  of  the  S0A_W  program  (switch  the 
positions  of  the  "Phi"  and  "Dummy"  variables) . 
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Additional  comments  on  the  responses  to  the  14  questions  of  Step 
7  follow: 

1.  This  is  the  start  time. 

2.  More  than  one  objective  analysis  will  produce  a  series  of 

forecasts . 

3.  This  is  the  interval  in  seconds  from  the  start  time  (1.)  to 

the  first  forecast  and  each  forecast  thereafter.  The 
example  used  is  for  a  single  nowcast.  For  a  forecast,  the 
analysis  will  occur  first  at  time  (1.  +.3.)  and  continue  for 
(2.) 

4.  and  5.  The  limits  in  4.  and  5.  are  for  the  full  data  set 
that  will  be  used  in  the  objective  analysis.  The  units  are 
meters  for  4.  and  seconds  for  5.  To  include  data  both 
before  and  after  the  start  time  increase  5.  In  this  case, 
data  is  sampled  every  2  minutes  or  120  seconds,  therefore  by 
setting  5.  to  120  the  available  data  goes  from  6  values  at 
14:30:00  to  18  values  at  14:28:00,  14:30:00,  and  14:32:00. 

6.  This  is  limit  in  meters  and  seconds  for  the  elimination  of 

the  distant  observations  from  each  calculation.  These 
values  should  be  smaller  than  4.  and  5. 

7 .  The  maximum  number  of  influential  points  used  for  each 

calculation  and  within  the  limits  of  6.  As  this  values 
becomes  large  more  outlining  data  point  are  used.  This 
greatly  increases  the  computer  time  required,  without  a 
significant  increase  in  the  quality  the  final  result.  With 
too  few  influential  points  (certainly  3  or  less)  '  learest 
values  are  only  used  and  information  from  other  close  points 
are  lost  to  the  calculation.  At  the  limit  of  1,  the  single 
closest  point  is  used,  resulting  in  a  field  that  is  a  series 
of  "plateaus". 

8.,  9.,  &  10.  The  correlation  function  used  here  falls  off  from 
1  at  zero  by  e  to  -R  power,  where 
R  =  sqrt [ (delta_x/xfold) 2  +  (delta_y/yfold) 2  + 

(delta_t/tfold) 2] . 
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11.  and  12.  The  two  input  files.  The  data  file  and  the  file  of 
the  interpolated  positions. 

13.  and  14.  The  two  output  files.  This  first  is  the  forecasted 
field  of  the  values  and  the  second  is  the  error  field.  The 
expected  error  field  is  a  measure  of  the  confidence  at  each 
interpolated  point.  The  closer  the  data  point  are  to  the 
interpolated  point  the  greater  the  confidence  and  smaller 
the  expected  error. 


Chapter  3 
TEST  RUNS 


3 . 1  INTRODUCTION 

Objective  analysis  has  been  successfully  applied  to  large  data 
sets  of  freely  drifting  oceanographic  buoys.  However,  what 
occurs  to  the  results  when  the  data  available  to  the  objective 
analysis  program  is  systematically  reduced  to  the  limit  of  a 
single  available  point?  Does  the  technique  fail  catastroph¬ 
ically,  or  do  the  results  simply  degrade  with  fewer  data  availa¬ 
ble  for  input?  And  if  the  results  just  degrade,  how  do  they 
degrade?  To  answer  these  questions,  a  series  of  test  runs  of 
the  scalar  objective  analysis  (SOA)  program  were  done.  The  data 
set  used  to  test  the  SOA  was  from  the  Fort  Pierce,  FL,  March- 
April  1985  Leeway  Experiment.  During  the  afternoon  of  1  April, 
six  surface  drifters  were  deployed  and  were  tracked  with  the 
Microwave  Tracking  System  (MTS) .  Lieutenant  Louis  Nash,  USCG, 
provided  the  edited  u  and  v  velocity  data  in  the  x-y  state 
coordinate  plane.  Figures  4-9.  Thus  we  had  a  data  set  which  was 
large  by  our  standards  for  investigating  the  SOA  program.  To 
further  simplify  the  problem  we  only  used  velocities  at  a  single 
point  in  time  -  14:30. 
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3 . 2  TEST  RUNS 


During  the  test  runs  of  the  objective  analysis  program  we  first 
tried  to  calculate  the  correlation  function  using  the  COM_COR 
program  on  all  six  drifters.  However,  six  drifters  proved  to  be 
insufficient  for  the  calculation  of  a  sensible  correlation 
function.  Therefore  we  used  the  analytic  form  of  an  exponential 
decay  function  as  given  by  equation  (10) ,  the  correlation 
function. 

Equation  (10)  was  used  as  the  basis  for  generating  current  field 
values  on  a  grid  of  15  by  15  regular  spaced  points  from  the 
observed  surface  current  data.  That  grid  has  points  between 
current  data  points  and  to  a  limited  extent  outside  the  region. 
The  question  was  then  asked,  "What  is  the  minimum  number  of 
drifters  needed  by  the  SOA  program  to  produce  a  sensible  current 
field?”.  To  start,  all  six  drifters  were  used  to  define  a 
"best"  estimate  of  the  existing  current  field.  This  provided 
the  best  estimate  of  the  u-current  field  against  which  the  other 
fields  were  compared.  When  drifters  were  removed  one  at  a  time 
from  the  data  set,  the  forecasted  fields  changed.  The  fields 
generated  by  the  reduced  number  of  drifters  were  then  compared 
to  the  results  with  the  "best  estimate"  field  determined  by  all 
six  drifters.  Two  different  schemes  for  removing  drifters  were 
investigated.  Case  "A"  was  to  remove  the  least  useful  drifters 
one  at  a  time  leaving  the  "best"  to  be  used  to  estimate  the  U- 
field.  Case  "B"  was  to  remove  the  most  useful  drifters  one  at  a 
time  leaving  the  "worst"  to  provide  the  estimate  of  the  U-field. 
Both  schemes  are  shown  in  Table  I. 
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TABLE  I  -  DRIFTER  CASES 


CASE-A  DRIFTERS  CASE-B  DRIFTERS 

(BEST) _ (USED) _ (WORST) _ (USED) 


6 

5, 

6,  7, 

8,  9, 

10 

6 

5, 

6, 

7, 

8,  9, 

10 

5A 

5, 

7, 

8,  9, 

10 

5B 

5, 

6, 

7, 

9, 

10 

4A 

7, 

8,  9, 

10 

4B 

5, 

6, 

7, 

9 

3A 

7, 

8, 

10 

3B 

6, 

8,  9 

2A 

8, 

10 

2B 

6, 

9 

1A 

8 

TABLE  I: 

The 

Best 

(A) 

and  the 

Worst 

(B) 

cases 

used 

for 

comparisons  of  subsets  of  drifters  with  the  full  set 
of  six  drifters. 


The  differences  were  calculated  at  each  15  by  15  grid  point 
between  the  forecasted  field  derived  from  all  six  drifters  and 
each  15  by  15  grid  point  of  the  cases  above.  The  differences 
from  all  the  225  points  produced  the  measured  error  fields.  The 
square  root  of  the  sum  of  the  squared  differences  is  an  index  of 
the  cumulative  measured  error.  The  equation  for  cumulative 
measured  error  index  for  case  5A  is  shown  below: 


Cumulative  Measured  2 

Error  Index  ^  i j  )  (11) 

where  the  subscripts  "ij "  are  the  elements  of  the  15  by  15 
array.  Therefore,  we  have  fields  of  the  forecasted  U-velocity, 
expected  error,  measured  error  and  an  index  of  the  cumulative 
error. 

Figure  10  shows  the  eastward  component  of  velocity  (u)  at  14:30 
from  the  six  drifters  and  Figure  11  is  the  northward  (v) 
component  of  velocity.  In  Figure  12  the  current  field  generated 
by  the  SOA  scheme  for  the  six  drifters  is  shown.  The  currents 


near  the  drifters  are  reproduced  faithfully  enough,  including 
the  region  of  currents  less  than  10  cm/s  in  the  upper  central 
region  and  the  currents  greater  than  16  cm/s  in  the  southwest 
corner.  Figure  13  shows  the  expected  errors  and  they  are  low 
across  the  region  with  slightly  higher  values  in  the  four 
corners,  as  expected.  The  measured  error  for  the  6  drifters 
case  is  zero  by  definition,  and  therefore  is  not  shown.  The 
objective  analysis  forecast  of  the  U-fields  for  scheme  "A"  are 
presented  in  Figures  14  -  17.  The  expected  error  generated  by 
the  objective  analysis  program  for  the  "A"  scheme  are  shown  in 
Figures  18  -  21.  The  fields  of  measured  error  for  scheme  "A" 
are  shown  in  Figures  22  -  25.  The  "B"  scheme  forecast  fields, 
expected  error  fields  and  measured  error  fields  are  presented  in 
Figures  26  -  37. 

The  forecast  fields  for  all  scheme  "A",  Figures  14  -  17, 
maintains  the  basic  field  forecasted  by  the  full  six  drifters, 
Figure  12.  This  is  also  shown  by  the  measured  error  fields, 
Figures  22  -  25,  which  are  generally  low  right  across  the  field. 
Even  "2 A"  which  uses  just  drifters  #8  and  #10  reproduces  the  U- 
field  without  trouble.  Figures  17  and  25.  This  is  in  direct 
contrast  with  scheme  "B"  forecast  fields,  Figures  26  -  29.  Here 
we  can  see  that  even  with  five  drifters  the  forecasted  field 
differs  greatly  from  the  original  field,  Figure  12.  The 
measured  error  fields  confirm  this,  Figures  34  -  37. 

The  cumulative  measured  error  index  for  the  two  schemes  are 
presented  in  Table  II. 
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TABLE  II. 

CUMULATIVE  MEASURED  ERROR  INDEX 


Number 

of  drifters  used  for 

Objective 

Analysis 

6 

5  4 

3 

2 

1 

Scheme  A 

0.0 

3.2  5.7 

13.3 

14.9 

49.4 

Scheme  B 

0.0 

14.7  18.9 

27.3 

24.2 

— 

TABLE  II.  The  Cumulative  Measured  Error  Index  for  Scheme  A 
and  Scheme  B. 


This  suggests  that  "2A"  did  as  good  a  job  as  predicting  the  U- 
field  as  did  "5B''. 

Scheme  "1A"  with  just  drifter  #8  is  the  special  case  of  the 
objective  analysis  taken  to  its  limit.  The  forecast  field. 
Figure  38,  is  uniform  at  the  value  of  the  drifter  (9.028  cm/s). 
The  expected  error  field,  Figure  39,  increases  away  from  the 
location  of  the  drifter  due  to  the  decrease  in  the  correlation 
function,  Equation  11.  The  measured  error  field.  Figure  40  has 
the  same  form  as  the  original  forecast  field  for  the  6  drifter 
(Figure  12) ,  but  offset  by  9  cm/s. 

3 . 3  RESULTS.  OF  TEST  RUNS 


The  main  portion  of  the  objective  analysis  technique  did  not 
fail  catastrophically.  Instead,  the  resulting  velocity  fields 
decreased  relative  to  the  original  field,  reading  from  left  to 
right  -  Table  II.  However,  the  more  important  point  is  that  the 
placement  of  the  buoys  or  drifters  relative  to  the  flow  field 
features  will  greatly  influence  how  well  the  total  flow  field  is 
estimated,  reading  from  top  to  bottom  -  Table  II.  An  example  of 
an  application  would  be  the  placing  of  two  buoys  on  either  side 
of  a  oceanic  front,  estimating  the  currents  better  than  several 
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buoys  all  on  the  same  side  of  the  front.  A  few  well  placed 
buoys  will  provide  better  information  than  many  poorly  placed 
buoys . 

The  determination  of  the  correlation  function  (Section  2.2)  from 
the  statistics  did  fail  for  the  data  set  used  here.  The  number 
of  buoys  required  to  determine  the  correlations  functions  is 
greater  than  six.  Therefore  an  analytical  formula  (Equation  10) 
was  used  instead.  In  this  study  neither  the  v-velocity  (north  - 
south)  field  nor  the  time  dependent  field  were  considered. 
Since  there  should  be  a  relationship  between  the  u  and  v  fields 
and  between  the  spacial  and  temporal  fields,  correlations  in 
time  could  be  used  to  infer  correlations  in  space.  Thus,  the 
buoys '  track  histories  could  be  used  to  determine  the  time 
correlation  which  then  would  provide  a  better  estimate  of  the 
velocity  field.  Further  work  is  being  conducted  on  improving 
estimates  of  the  current  field  from  limited  amounts  of  drifter 
data. 

CHAPTER  4 
CONCLUSIONS 

4 . 1  CONCLUSIONS 

The  search  planner  requires  real-time  surface  current 
information  to  determine  the  most  probable  location  of  the 
survivors  and  survivor  craft.  The  displacement  with  time  of  the 
search  datum  can  be  estimated  from  the  track  histories  of 
freely-drifting  surface  buoys.  All  types  of  freely  drifting 
buoys  will  generate  data  sets  that  consist  of  a  time  series  of 
irregular-spaced  positions.  As  a  set  of  buoys  freely  drift  in 
the  ocean,  (1)  positions  at  any  one  time  are  not  on  a  regular 
grid,  (2)  the  displacement  from  the  previous  positions  of  that 
buoy  to  the  next  position  is  uneven,  and  (3)  the  relative 
displacement  among  all  the  buoys  is  changing  with  time. 
Additionally,  some  types  of  buoys  (e.g.  TIROS  satellite  track 
buoys)  positions  are  reported  intermittently.  The  results  is  a 
"messy”  data  set. 
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The  results  of  the  test  runs  suggest  that  a  few  well  placed 
buoys  can  provide  a  better  estimation  of  the  current  field  than 
many  poorly  placed  buoys.  Therefore,  the  use  of  remote-sensing 
information  (e.g..  Side  Looking  Airborne  Radar  (SLAR)  derived 
sea  surface  roughness  or  N0AA6  satellite  Advanced  Very  High 
Resolution  Radiometer  (AVHRR)  derived  sea  surface  temperature 
imagery)  will  aid  in  the  determination  of  the  optimal  placement 
of  buoys. 

To  take  advantage  of  computer  analysis  of  the  buoys  tracks,  the 
data  set  must  first  be  cleaned  up.  The  essence  of  objective 
analysis  is  that,  given  an  irregular-spaced  data  set  in  time  and 
space,  a  data  set  is  interpolated  and  extrapolated  to  a  regular 
(even  and  square)  grid.  From  this  many  standard  computer 
programs  are  available  that  can  be  brought  to  bear  on  the 
analysis  of  the  buoy  tracks.  This  analysis  can  then  provide  the 
most  information  in  a  timely  manner  for  the  search  planner. 

In  Figure  41  the  Search  and  Rescue  Problem  is  blocked  out  (A) 
and  the  Datum  Movement  components  of  SAR  Planning  are  expanded 
in  (B)  .  As  the  available  Environmental  Data  sets  grow  in 
number,  size  and  complexity,  the  Coast  Guard  will  have  to 
provide  a  system  to  handle  these  data  sets.  Figure  41(B) 
illustrates  the  role  Objective  Analysis  would  play  in  a  possible 
future  system.  Environmental  Data  in  the  form  of  winds  (e.g., 
FNOC,  remotely  reporting  winds  from  ships,  meteorological-buoys, 
weather  stations,  and  Next  Generation  RADAR)  and  surface 
currents  (e.g.,  Loran-C  and  ARGOS  buoys,  and  satellite  imagery) 
would  have  to  be  collected  and  processed  initially  before 
passing  to  the  Objective  Analysis  module.  Then  the  Objective 
Analysis  module  would  assimilate  and  convert  these  irregular 
spaced  data  sets  into  standard  arrays.  The  standard  arrays  of 
wind  would  provide  nowcasts  of  the  Wind  Field  for  the  Leeway 
Calculation.  The  standard  arrays  of  the  surface  currents  could 
either  be  the  nowcast  of  the  Ocean  Current  Fields  or  the  Initial 
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Conditions  required  by  the  numerical  model  of  the  local  physical 
oceanography.  The  numerical  model  would  then  provide  forecasts 
for  the  Ocean  Current  Field.  The  Ocean  Current  Field  along  with 
the  Leeway  Calculation  would  be  used  for  Datum  Movement 
calculations . 

In  summation,  an  objective  analysis  program  was  modified  and 
tested  to  demonstrate  its  capabilities  on  small  data  sets.  The 
results  of  the  test  show  that  objective  analysis  can  effectively 
work  on  as  little  as  two  buoy  drift  tracks.  The  Coast  Guard  is 
going  to  utilize  a  variety  of  environmental  data  sources  to 
effectively  locate  the  most  probable  positions  of  survivors. 
The  objective  analysis  module  will  be  a  first  step  in  preparing 
these  data  sets  for  use  by  the  search  planner.  The  application 
of  these  data  sets  to  other  computer  programs  will  provide  the 
search  planner  with  the  most  probable  positions  of  the 
survivors.  Work  is  continuing  on  the  objective  analysis  module. 

4.2  FUTURE  WORK 

Presently  the  entire  front  end  and  back  end  of  this  set  of 
programs  contain  room  for  improvements.  Work  has  started  to 
greatly  increase  the  ease  of  inputting  diverse  data  sets, 
including  data  sets  where  positions  are  given  in  latitude  and 
longitude.  Further  work  is  planned  to  use  objective  analysis 
for  a  variety  of  data  sets:  (1)  the  wind  velocity  fields  using 
input  data  from  several  surrounding  weather  stations,  (2) 
surface  current  fields  from  NOS  tidal  data,  (3)  results  of  local 
numerical  models,  (4)  satellite  AVHRR  sea-surface  temperature 
imagery  data,  (5)  sea  surface  currents  determined  from 
successive  AVHRR  imageries,  and  (6)  the  relatively  large  buoy 
data  set  collected  by  the  Commander,  International  Ice  Patrol 
( CUP)  .  The  CUP  data  set  will  provide  experience  with  a  data 
set  similar  to  what  a  Group  or  District  would  collect  after 
several  years  of  deploying  VHF  Loran-C  or  TIROS  type  datum 
marker  buoys. 
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FIGURES 


Figure  1.  The  defined  correlation  function  (equation  10),  where 
z  is  the  value  of  the  function  and  xfold  and  yfold  -1.  The 
circle  at  z  =  0.368  is  the  e-folding  scale. 


Figure  2.  The  defined  correlation  function  (equation  10  )  where 
z  is  the  value  of  the  correlation  function  and  xfold  =  1  and 
yfold  =  3.  The  ellipse  at  z  =  0.368  is  the  e-folding  scale. 
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DRIFTER  10 


Figure  9.  The  current  speeds  of  Drifter  No.  10  as  a  function  of 
time.  The  solid  and  dotted  lines  indicate  the  u-  and  v- 
components  of  current  speeds  respectively. 
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Figure  10.  Position  and  u-component  speeds  at  14:30,  1  April 
1985.  The  symbols  "+M  indicate  the  drifter  location.  The 
numerics  above  the  M+"  denote  the  drifter  ID;  while  the  numerics 
below  "+"  denote  the  values  of  the  u-component  (cm/sec) . 
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Drifter  Position  &  V_speed 
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Y  of  f  : 349000 
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Figure  11.  Position  and  v-component  speeds  at  14:30,  1  April 
1985.  The  symbols  "+"  indicate  the  drifter  location.  The 
numerics  above  the  "+"  denote  the  drifter  ID;  while  the  numerics 
below  "+"  denote  the  values  of  the  v-component  (cm/sec) . 


0.  fl.  Forecast  Field 
Time:  1  Rpr  1985  14:30:00 


Figure  12.  Objective  Analysis  Forecast  Field  of  the  u-component 
of  velocity  for  Case  6. 
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Min.:  .034;  Max.:  .243;  Cl:  .02 
DRIFTERS  5,  7,  8,  9,  10 

Figure  18.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  5A. 
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Expected  Error  Field 
Time:  1  Rpr  1985  14:30:00 


Figure  19.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  4A. 


Expected  Error  Field 
Time:  1  Rpr  1985  14:30:00 


Figure  20.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  3A. 
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0.  fl.  Forecast  Field 
Time:  1  flpr  1985  14:30:00 
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0.  R.  Forecast  Field 
Time:  1  Rpr  1985  14:30:00 


Figure  29.  Objective  Analysis  Forecast  Field  of  the  u-component 
of  velocity  for  Case  2B. 


Expected  Error  Field 
Time:  1  Rpr  1985  14:30:00 


Figure  30.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  5B. 
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Expected  Error  Field 
Time:  1  Apr  1985  14:30:00 


Figure  31.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  4B. 


Expected  Error  Field 
Time:  1  Apr  1985  14:30:00 


Figure  32.  Expected  Error  Field  of  the  u-component  of  velocity 
for  Case  3B. 
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Figure  34.  Measured  Error  Field  of  the  u-coxnponent  of  velocity 
for  Case  5B. 
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Min.:  -.23}  M«x.:  3.129;  Cl:  .5 

MEASURED  ERROR  of  (6  -  3B)  I*  27.3 


Figure  36.  Measured  Error  Field  of  the  u-component  of  velocity 
for  Case  3B. 
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of  velocity  for  Case  1A,  Drifter  No.  8, 
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SEARCH  AND  RESCUE  PROBLEM  DEFINITION 


(A) 


DATUM  MOVEMENT  FOR  SAR  PLANNING 


Figure  41.  The  Search  and  Rescue  Problem  Definition  (A) 
and  the  Datum  Movement  for  Search  and  Rescue  Planning  (B) 
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APPENDIX  A 

DEFINITIONS  OF  PROGRAM  VARIABLES 

Three-dimensional  array  of  the  correlation 
function. 

Observed  data  file. 

Date  and  Time  of  the  analysis  to  make. 

Time  interval  for  each  time  extrapolation. 

Time  lag. 

X-  and  y-spatial  lags. 

The  error  variance  of  the  u-component  current 
speed. 

The  error  variance  of  the  v-component  current 
speed. 

Expected  error. 

Error  variance. 

Maximum  number  of  influential  points. 

Center  of  the  correlation  table. 

Maximum  space  lags. 

Maximum  time  lags. 

Number  of  input  observed  data. 

Number  of  inter/extrapolation  positions. 

Number  of  objective  analysis  to  make. 

Number  of  data  used  in  the  analysis. 

Optimal  data  used  in  the  analysis. 

Restricted  observational  data. 

Expected  error  output  file  for  scalar  objective 
analysis. 

Output  file  for  scalar  objective  analysis 
forecast  fields. 

Time  interval  in  the  correlation  look-up  table. 
Forecast  time. 

Estimated  value. 

Time  limit  before  the  forecast  time. 

Maximum  time  radius  before  and  after  the  time 
of  the  analysis. 

Observation  time  of  optimal  data. 

Observational  time  of  restricted  data. 

Time  limit  after  the  forecast  time. 

The  correlation  function  between  the  current 
components . 

The  error  variance  of  u-  and  v-component 
current  speeds. 

The  estimated  value  of  u-  and  v-component 
current  speeds. 

Expected  error  output  file  for  vector  objective 
analysis. 

Output  file  for  vector  objective  analysis 
forecast  fields. 

X-  and  y-space  intervals  in  the  correlation 
look-up  table. 

X-  and  Y-coordinates  of  interpolation 
positions. 
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Xlimit 

Xopd, Yopd 
Xrd, Yrd 


APPENDIX  A  (cont'd) 

DEFINITIONS  OF  PROGRAM  VARIABLES 

Maximum  spatial  radius  from  the  reference  point 
of  domain. 

X-  and  Y-coordinates  of  optimal  data. 

X-  and  Y-coordinates  of  restricted  data. 


A-2 


APPENDIX  B 


PROGRAM  LISTINGS 

Page 

PROGRAM  COM_COR .  B-3 

SUBPROGRAM  FNCov_f .  B-6 

SUBPROGRAM  FNCOV_t .  B-7 

PROGRAM  GET_DATA .  B-8 

PROGRAM  LINDAT_UV .  B-9 

PROGRAM  PLOT_PUS .  B-13 

PROGRAM  PLOT_PVS .  B-17 

PROGRAM  PREPARE .  B-21 

PROGRAM  PRINT_OA .  B-22 

PROGRAM  GET_IP .  B-24 

PROGRAM  SOA_W .  B-25 

PROGRAM  VOA_W .  B-43 

PROGRAM  SOA_DIFF .  B-61 

PROGRAM  PLOT  OA .  B-6 2 
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[  BLANK  ] 
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1000  I  Program  COM_COR 
1005  I 

1010  !  -  routine  computes  the  correlation  of  observed  data 

1015  !  as  a  function  of: 

1020  !  Delta_x  and  Qelta_y  (space  lags) 

1025  !  Delta_t  (time  lag). 

1030  ! 

1035  OPTION  BOSE  1 

1040  DIM  Xdata( 1089 ) ,Ydata( 1089 >. Tdata<  1089 > ,Phi( 1089 > 

1045  DIM  Bina( 33,33) ,Binc( 33,33) ,Binx< 33,33) ,Biny( 33,33) 

1050  ! 

1055  Clear$=CHR$(  255 )&"K“ 

1060  OUTPUT  2  USING  "#,K" ; Clear®  !  clear  the  screen. 

1065  PRINT  TABXY( 10, 10) ."Enter  the  number  of  bins  in  x -  and  y-coordinates . " 
1070  INPUT  "Npx? ,Npy?" ,Npx ,Npy 

1075  OUTPUT  2  USING  "# , K " i Cl  ear*  !  clear  the  screen. 

1080  PRINT  TABXY( 10, 10) ."Enter  x-  and  y-bin  sizes  and  time  lag:" 

1085  INPUT  "X_bin? ,  Y_bin?,  T_bin?" , X_bin , Y_bin , T_bin 
1090  OUTPUT  2  USING  , K " ; Cl  ear*  !  clear  the  screen. 

1095  PRINT  TABXY<  10, 10) ."Enter  number  of  iterations" 

1100  INPUT  " I terat ions?" , Interat ions 

1105  OUTPUT  2  USING  "#,K" ; Clear®  !  clear  the  screen. 

1110  PRINT  TABXY( 10, 10) , "Enter  the  correlation  filename" 

1115  INPUT  "Corr_f ile$?“ ,Corr_f ile$ 

1120  OUTPUT  2  USING  "#,K";Clear$  !  clear  the  screen. 

1125  PRINT  TOBXY(  10, 10)  .“Enter  the  observed  data  file'¬ 
ll  30  INPUT  "Data_f ile$?" ,Data_f ile® 

1135  OUTPUT  2  USING  "#,K" ;Clear$  !  clear  the  screen. 

1140  ! 

1145  1 - get  data 

1150  ! 

1 155  COLL  Get_data<  Data_f i 1 e$ , Xdata( *  > ,Ydata<  * ) ,Tdata(  * ) ,Phi( * ) .Ndata) 

1 1 60  Mean=0. 

1165  FOR  1=1  TO  Ndata 
1170  Mean=Mean+Phi(  I ) 

1 175  NEXT  I 

1180  Mean=Mean/Ndata 

1185  PRINT  "Number  and  mean  of  input  data:  ".Ndata, Mean 
1190  i 

1195  !  -  remove  the  mean  of  input  data 

1 200  ! 

1205  FOR  1=1  TO  Ndata 

1210  Phi< I )=Phi< I )-Mean 

1215  NEXT  I 

1220  I 

1225  Irep=0 

1230  T_inc=0. 

1235  M_bin=<  Npx  +  1  )/2 
1240  N_bin=( Npy+ 1  )/2 
1245  REPEAT 
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1Z50  T_inc=T_i nc+Tb 1 n* I  rep 

1255  Irep=Irep+l 
1250  FOR  1=1  TO  Npx 
1255  FOR  J=1  TO  Npy 

1270  Bina< I , J)=0. 

1 Z75  Bi nc( I , J ) =0. 

1280  Binx< I , J )=0. 

1285  Bi ny( I , J ) =0. 

1290  NEXT  J 

1295  NEXT  I 

1300  i 

1305  !  -  loop  to  find  pairs  of  data  separated  by  T_bin+-5  time  units. 

1310  ! 

1315  FOR  J=1  TO  Ndata 

1320  FOR  1=1  TO  Ndata 

1325  Delta_t  =  Tdata(  J)-Tdata(  I ) 

1330  IF  Delta_t>T_inc+ . 5*T_bin  THEN  GOTO  Next_i 

1335  IF  Oelta_t<T_inc- . 5*T_bin  THEN  GOTO  Next_i 

1340  Delta_x  =  Xdata(  J>~Xdata< I ) 

1345  Del ta_y=Ydata<  J  >- Ydata(  I ) 

1350  Idx=INT( ( Del ta_x+ . 5*X_bin ) /X_bin)+M_bin 

1355  Idy=INT<  (  Del ta_y+ . 5* Y_bin  > / Y_bin  > +N_bin 

1350  IF  Idx<1  OR  Idx>Npx  THEN  GOTO  Next_i 

1355  IF  Idy< 1  OR  Idy>Npy  THEN  GOTO  Next_i 

1 370  Binc< Idx , Idy )=Binc( Idx , Idy )+1 

1  375  Bina( Idx , Idy )=Phi(  J )*Phi< I )+Bina< Idx , Idy ) 

1380  Binx< Idx , Idy )=Phi< J )*Phi< J )+Binx( Idx , Idy ) 

1385  Biny( Idx , Idy )=Phi(  I )*Phi< I )+Biny< Idx, Idy) 

1390  Next_i : NEXT  I 
1395  NEXT  J 
1400  ! 

1405  !  -  normalize  bin 

1410  ! 

1415  BEEP 

1420  FOR  J=1  TO  Npy 

1425  FOR  1=1  TO  Npx 

1430  IF  Binx< I , J )=0.  OR  Biny(I,J>=0.  THEN  GOTO  Nexti 
1435  Bina( I , J)=Bina( I , J)/SQR(Binx( I , J)*8iny< I , J) ) 

1440  Nexti:  NEXT  I 
1445  NEXT  J 
1450  i 

1455  1  -  output  correlation  matrix 

1450  ! 

1455  ASSIGN  @Path_out  TO  Corr_file$ 

1470  FOR  1=1  TO  Npx 

1475  FOR  J=1  TO  Npy 

1480  OUTPUT  @Path_out;Bina( I ,J) 

1485  NEXT  J 

1490  NEXT  I 

1495  UNTIL  Irep=Interat ions 
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1500  Finish:  i 

1505  ASSIGN  @Path  out  TO  * 

1510  END 
1515  i 
1 5Z0  i 
1 5Z5  i 

1 530  Get_dat a: SUB  Get  _dat a<  F  i I enanel , Xdat a<  * ) , Ydata<  * ) ,  Tdata< * ) , Phi datai *  > , Ndata 
) 

1535  ASSIGN  @Path._in  TO  Filename® 

1540  Read_in_data:  L 

1545  Ndata=0 

1550  PRINT  "  **********  Echo  check  of  the  first  ten  records  ********** 

*  " 

1555  PRINT  "  No.  X_POS  Y_POS  TIME  DATA" 

15G0  ON  END  @Path_in  GOTO  1610 

1565  DISP  USING  "K ,K ,K" ; "Reading  the  observed  data  from  ".Filename®.",  please 
wait .  " 

1570  ENTER  @Path_in i X . Y . T . Phi , Dummy 

1575  IF  Ndata<= 10  THEN  PRINT  USING  "DDDD, XX , MD. DDDDDDE , XX ,MD. DDDDDDE , XX ,MD. DDD 

DDDDDDE . XX . MD . DDDDDDE " ;Ndata,X,Y,T,Phi 

1580  Ndata=Ndata+ 1 

1585  Xdata(Ndata)=X 

1530  Ydata< Ndata )=Y 

1595  Tdata( Ndata)=T 

1600  Phidata<Ndata)=Phi 

1605  GOTO  1560 

1610  ASSIGN  @Path_in  TO  »!  Closing  input__file. 

1615  DISP 
1 6Z0  SUBEND 
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i  000  i 

1005  !  SUBPROGRAM  FNCov_f 

1010  I 

1015  DEF  FNCov(Delx,Dely,Delt,Xfold,Yfold,Tfold) 

1020  i 

1025  !  routine  gets  the  correlation  function  from  a  fitted  formula. 

1030  ' 

1035  OPTION  BASE  1 
1040  RAD 

1045  COM  /Efield/  Evar 

1046  RETURN  EXP(  -SQR<  (  Del  x / Xf  ol d  )  *  (  Del x/Xf  ol  d )  +  (  Del  y/  Yf  ol  d  )  » <  Del  y/ Yf  ol  d )  +(  Del  t 
/Tfold>*< Delt/Tfold) ) ) 

1050  '  R2  =  SQR(  (  DelxAZ+Dely''Z  )/5.0E  +  4AZ  ) 

1055  i  Cor=EXP< -R2 ) 

1060  i  RETURN  Cor 
1065  FNEND 


B-6 


1000  i  SUBPROGRAM  FNCov_t 
1005  i 

1010  DEF  FNCov(  Del  ta  x , Del  ta_y , Del ta_t ) 

1015  i 

1020  1  -  routine  determines  the  correlation  function  from  the  look-up  table, 

1025  1  given  space  and  time  lags. 

1030  1 

1035  1  Corr  is  a  three-dimensional  array  of  the  correlation  function. 

1040  !  Delta_x  and  Delta_y  are  the  x-  and  y-spatial  lags. 

1045  !  Delta_t  is  the  time  lag. 

1050  1  M_bin  is  the  mid-point  in  the  x-dimension  of  correlation  table. 

1055  1  N_bin  is  the  mid-point  in  the  y-dimension  of  correlation  table. 

1060  1  T_bin  is  the  time  interval  in  the  correlation  look-up  table. 

1065  1  X_bin  is  the  x-distance  interval  in  the  correlation  look-up  table. 

1070  !  Y_bin  is  the  y-distance  interval  in  the  correlation  look-up  table. 

1075  OPTION  BASE  1 
1080  RAD 

1085  COM  /Co rr/  Corr ( 33 , 33 , 1 0 ) , X_bin , Y_bin , T_bin , M_bi n , N_bin 
1090  Idx  =  INT( ( Del ta_x+ . 5*X_bin ) /X_bin  >+M_bin 
1095  Idy=INT  < ( Del ta_y+ . 5* Y_bin ) / Y_bin > +N_bin 
1  100  Idt  =  INT(  (  Delta_t  +  . 5*T_bin  >/T_bin)  +  1 
1105  Cor=Corr( Idx , Idy , Idt ) 

1110  RETURN  Cor 
1115  FNEND 
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1000  i PROGRAM  GET_D ATA 
1005  i 

1010  DIM  X( 1089 ) ,Y( 1089 ), T<  1089 ) ,U<  1089 ) ,V< 1089 ) 

1015  i  DATA  “DRIFTER_5: ,700, 1 " ,0 
1020  '  DATA  "DRIFTER_G: ,700, 1" .0 

1025  DATA  "DRIFTER_7:  ,700, 1“  ,0 
1030  DATA  ”DRIFTER_8: .700. 1" .0 
1035  DATA  "DRIFTER_9: ,700 , I " ,0 

1040  DATA  "DRIFTER_10: ,700, 1" ,0 

1041  ON  ERROR  GOTO  I04G 

1042  CREATE  BDAT  “ OBSDAT4A: , 700 , 1" , 1 00 

1043  CREATE  BDAT  " S0A_FCST4A : , 700 , 1" , 1 0 

1044  CREATE  BDAT  " SOA_EVAR4A : , 700 , 1" , 1 0 
I 04B  OFF  ERROR 

1047  ASSIGN  @Pathout  TO  "0BSDAT4A: ,700, 1" 

1050  FOR  N= 1  TO  4 

1055  READ  F i lenameS , Nsk ip 

1060  PRINT  Filename® 

1065  Get_data<  F i 1 enameS , Nsk ip , X< * ) , Y<  » ) , T( # ) , U( * ) , V<  » ) .Ndata) 

1070  PRINT  Ndata 
1075  FOR  1=1  TO  Ndata 

1 080  OUTPUT  ePathout ;  X(  I  )  ,  Y<  I  )  ,  T<  I ) ,  U<  I  )  ,  V<  I  ) 

1085  NEXT  I 
1080  NEXT  N 
1095  END 

1 1 00  Get_data: SUB  Get_data(  F i lenameS ,Nsk ip , Xdata( * ) , Ydata( *  > , Tdata<  *  > ,  Udata( *  > , V 
data<  * ) , Ndat  : ) 

1105  DIM  Drno$l21 ,Id$l 1 1 ,Desc$l 301 
1110  INTEGER  Npos 
IMS  ASSIGN  @Path_in  TO  Filenames 
1120  ENTER  @Pat h_in; DrnoS , IdS , T$ , Npos 
1125  Ndata=0 

1130  ON  END  @Path_in  GOTO  1180 

1135  ENTER  @Path_in; X , Y , T  ,  U  ,  V 

1140  IF  Ndata< 1 0  THEN  PRINT  X,Y,T,U 

1145  Ndata=Ndata+ I 

1150  Xdata(Ndata)=X 

1155  Ydata< Ndata)=Y 

1160  Tdata(  Ndata )  =  T 

1165  Udata( Ndata ) =U 

1170  Vdata( Ndata >=V 

1  175  GOTO  1 1  30 

1180  ASSIGN  @Path_in  TO  *1  Closing  input_file. 

1185  SUBEND 
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I  000  >  Pr  01)1  .lit  I  I  NIIHI  UV 

1 00!.  I 

1010  1  Itinil  i  mi!  | j  1  < > t .  t  he  ii  .mil  v  cnmpoiicn t  current  speeds  ns  a 

I0l'i  1  f  urn :  l  i  on  of  tme.  I  he  solid  mill  Hotted  line:,  indicate  the 

10/0  1  ii  mid  v  component  speeds  respectively. 

10/G  i 

10  w>  op;  ion  until  i 

10.11,  DIM  Oil,  it.  i(  10FIII  )  .  Vil.it, i(  I0fl'l  ).  I(  I0HII  ) 

1040  hi  LOCH  It  H$l  I  0  I 
I  04 G  1  I  i  t  1  e  of  tins  plot 
I  0‘i0  IIHIH  "DM  1 1  ICR  ')" 

!0L,r,  l  X  axis  Label 
I0H0  IIHIH  "1  HPH  IHfl'," 
l0hS  1  Y  ii/s  i  s  l  ptie  1 
1070  IIHIH  "cm/ sec" 

I07D  l  X  Ilfl.X  RIliHI.X  tick  spar  i  nq  ,  X  pm  j  or  count 

1 0H0  OHIO  I,  /4I.  G 

l0H‘i  i  Y  bottom.  Y  top.Y  tick  sp, icing, Y  ma.|  or  count 
10110  IIHIH  0.  .  /ct.  ,  I  ,  !i 

I0!P>  •  Dump  graph l  c::>  to  printer  (yes/rio;  1/0) 

1100  DHIH  0 

I  I  0h  •  Dri  t  ii  s  timp  ling  i  ri  l  er  v.i  t  (  seconds  ) 

1110  DHIH  170. 

IIID  1  the  origin  of  time  axis 

I  1/0  DHIH  "  I  HP  Ft  I  98Fj“ 
ll/!i  DHIH  "10:00:00" 

1130  •  input  file  name,  number  of  data  to  be  skipped 

11.31,  DHIH  "DRIFTER  1. :  . 700 ,  I  "  , 0 
1140  C$T:HR$<  ?55>lt"K" 

11  4S  RE  HD  1 1 1 1  e$ 

I  I S0  RE  HD  X  1 abel* 

Mb'.,  RI:  HD  Y  label* 

1  ID0  ON  KBD  GOTO  Exit 

I  I  EiG  REHD  X  left.X  right, X  tick  spacing, X  major  count 
1170  REHD  Y  bottom, Y  top.Y  tick  spacing, Y  major  count 
117',  REHD  Dump  graphics 
I  I  H0  Major  tick  size" 3.0 

I  I  Fib  OUFPtir  /  USING  "tt.K'TC*  '  Clear  screen  for  graph 

I  1110  RF  HD  Del  ta  t 
I  I  Sb  RF  HD  Uat  e0  * 

1/00  RFHD  I i me®  * 

I  / 0G  I ime0»UHLE< Date0  * )  +  \  I  ME<  l ime0  $ ) 

1/10  RIHI)  Input  fileS.Nskip 

l/IF-i  Get  datat  l  npu  t  f  i.  L  e*  ,  Nsk  i  p  ,  Udatat  *  )  ,  Vdatat  »  )  ,  K  »  )  ,  Ndata  )  1  Get  data  for  p 

lotting 


7  7  V) 

//'j 

OlJIPUf  7 

gin i  r 

'  USING 

"tt.K'TC* 

1  Clear  screen  for  graph 

1  Initialize  various  graphics  parameters 

/  30 

iri  or  h  r 

IS  3," 

INI!  RNfiL" 

1  Use  the  internal 

screen 

/  3 1 

ess 

PL  Of  IF  ft 

70b 

IS  70! i 

, " HPGL  " 

1  Use  the  external 

HPGL  plotter  at  HP  IB 
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I Z 35  GRAPHICS  ON  1  Turn  on  the  graphics  screen 

1Z40  LORG  S  1  Reference  point:  center  of  top  of  label 

1245  1  Determine  hou  many  GDUs  wide  and  high  the  screen  is 

1  250  6du( Xgdu  max . Y_gdu_max ) 

1Z55  FOR  I=-.Z  TO  .2  STEP  .1  1  Offset  of  X  from  starting  point 

1ZG0  MOVE  X_gdu_max/Z+I , Y_gdu_max  1  Move  to  about  middle  of  top  of  screen 

1ZG5  LABEL  USING  j Title®  1  Write  title  of  plot 

1 Z70  NEXT  I  1  Next  position  for  title 

1275  DEG  !  Angular  mode  is  degrees  (used  in  LDIR) 

IZ80  Label  (  4 . 8  ,  .  6 , 90,  S  ,  I  ,  0.  ,  Y_gdu_.max/Z  ,  Y_1  abel$  )  !  Write  Y_axis  label 

1Z85  Label  (  3 . 8  ,  .  6 ,0 , 4 , 1  ,  X_gdu_max  /  Z  ,  .  01  *  Y_gdu_max  .  X_1  abel  $  )  1  Write  Xaxis  label 


1Z90  1  Define  subset  of  screen  area 

1  Z95  VIEWPORT  . 1 *X_gdu_max , . 9*X_gdu_max , . 1  * Y_gdu_max , . 9  * Y_gdu_max 

1300  1  Anisotropic  scaling:  lef t/nght/bottom/top 

1305  WINDOW  X_1 eft , X  r lght , Y_bot t om , Y_t op* 1 .01 

1310  1  Draw  axes  intersecting  at  lower  left 

1315  Y_ax is_locat ion-X_lef t 

13Z0  X_axis_locat ion=Y_botton 

1 3Z5  AXES  X_t ick_spac ing , Y_t ick_s pacing , Y_ax is_l ocat ion , X_ax is_l ocat ion , X_maj  or 
_count , Ymaj  or  count , Maj  or_t ick_size 


1  330 
1335 
1  340 
1  345 
1  350 
1  355 
1  3G0 
1  3G5 
I  370 
1  375 
1  380 
1  385 
1  390 
1  395 
1400 
1405 
1410 
14)5 
14Z0 
1  4Z5 
1430 
1435 
1440 
1445 
1450 
1455 
1460 
1  4G5 
1470 
1475 


CLIP  OFF 
CSIZE  3.5, .6 
LORG  G 

WINDOW  X_left ,X_right , . 


So  labels  can  be  outside  VIEWPORT  limits 
Smaller  chars  for  axis  labelling 
Ref.  pt :  Top  center  !\ 

»  Y_gdu_max , . 9  *  Y_gdu_max 


X_step=X_t ick_spacing*X_maj  or_count 
FOR  I  =X_1  ef t  TO  X_nght  STEP  X_step  1 
MOVE  I-X_left+1  , . 09* Y_gdu_max I  A  smidgeon  below  X-axis 
A$=TIME$(Time0+< 1-1 )*Delta_t> 

LABEL  USING  "tt,K"5A$l 1 ; 2 ] 

NEXT  I 

WINDOW  X_1 ef t , X_r ight , Y_bottom , Y_top* 1 .01 

LORG  8  !  Ref.  pt :  Right  center 

Y_step=Y_t ick_spacing* Y_maj  or_count 

FOR  I=Y_bottom  TO  Y_top  STEP  Y_step  ! 


\ 


>  Label  X-axis 


MOVE  .8,1 

LABEL  USING  "#,K"sI 
NEXT  I 
PENUP 

!  Anisotropic  scaling: 


Smidgeon  left  of  Y-axis 
DD.D;  no  CR/LF  !  / 

et  sequens 

LABEL  statement  leaves  the  pen  down 
lef t /right /bottom/ top 


\ 


>  Label  Y-axis 


WINDOW  X_lef t ,X_right , Y_bottom , Y_top* 1 .01 
Is=INT( T (  1  )-Time0)/Delt  a_t  + 1 
LINE  TYPE  1 

FOR  1=1  TO  Ndata  !  Points  to  be  plotted... 

PLOT  I +Is , Udat a( I )  I  Get  a  data  point  and  plot  it  against  X 

NEXT  I 
PENUP 

LINE  TYPE  8 

FOR  1=1  TO  Ndata  1  Points  to  be  plotted... 

PLOT  I+Is , Vdata< I )  1  Get  a  data  point  and  plot  it  against  X 

NEXT  I 
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I4B0  1  T)[SP  "Inter  'Space  liar  ’  to  qo  ori" 

I486  UP  I  I  I 
1490  DISP 

1495  GOTO  148(5  1  View  the  plot  ns  long  as  you  uant 

I 500  1  a  1 1  ■  If  Dump  graph i cs  I HEN 

1606  '  DUMP  GRAPHICS  CRT  10  #701 

IS  10  END  IF 

IMG  GRAPHICS  OIF 

I 5Z0  OUTPUT  a'  USING  "tt,K";C$ 

I 5Z5  FND  i  finis 

1530  i 

1535  Gdu: SUB  Gdu(  X  gdu  max  ,  Y  qdu  nax  ,  OPTIONAL  Gdu  x  pi  i  d  ,  Gdu  ymid) 

I 540  1  This  returns  Kright.  Yhiqh  and  their  respective  midpoints  in  GDUs. 

1545  1  Note  that  if  Gudxmid  is  defined,  Gdu  ymid  must  be  also. 

1550  COM  /Gunits/  Gduxmax , Gduymax , Udu  xmi n , Udu  xmax , Udu  ymi n , Udu  ymax , Show 
1 555  IF  Gdu_xmax«0  THEN 
1 560  Gdu_xmax=100»MAX< I .RATIO) 

1565  Gdu _y max  =  1 00* MAX (  I  , I /RA  F10 ) 

1570  END  IF 

1575  X  gdumax =Gdu_xmax 

1580  Y  gdu_max=6du_ymax 

1585  IF  NPAR>2  THEN 

1580  Gud  xmid=Gdu. xmax * . 5 

1585  Gud..ymid’"Gdu_ymax  »  .  5 

I  G00  END  IF 

IG05  SUBEND 

1GI0  Label: SUB  Label t Cs lze ,Asp_rat io , Ldir , Lorq .Pen , X , Y ,  lex t$ ) 

1615  1  This  defines  several  systems  variables  (in  CSIZE,  LOIR,  etc.),  arid 

I6Z0  1  labels  the  test  (if  any)  accordingly. 

IGZ5  DEG 

1630  CSIZE  Cs  tze  ,  Asp__rat  io 

1635  LDIR  Ldir 

1640  LORG  Lorg 

1645  PEN  Pen 

1650  MOVE  X  ,  Y 

1655  IF  TextSO""  THEN  LABEL  USING  "#,K";Text$ 

1660  PENUP 

1665  SUBEND 
1670  i 

1  675  Ge  t._data:  SUB  Get__data<  File_in$,Nskip,Udata(»),  Vdata(  *  ) ,  T<  *  )  ,  Ndata) 

1680  I 

1685  OPTION  BASE  1 

1690  DIM  DrnoSI Z 1 , I d$I 11 , Oesc$I 301 

1695  INTEGER  Npos 

1700  ASSIGN  @Path_ l n  TO  File_in* 

1705  'read  in  the  header 

1710  ENTER  @Path_ in; Drno* , I d$ , Desc* , Npos 
1715  IF  Nsk  ip=0  THEN  GOTO  Read._i n _dat a 
1 7Z0  FOR  N“ I  TO  Nsk ip 
I  7Z5  ENTER  @Path_ in ; Xd , Yd , Td , Ud . Vd 
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1730  NEXT  N 

1735  Read_in_data: 1=  I 

1740  ON  END  @Path_in  GOTO  1760 

1745  ENTER  @Path_ in ; X , Y , T< I > , Udata< I ) , Vdata< I > 

1750  1=1+1 

1755  GOTO  1740 

1750  Ndata=I-1 

1765  ASSIGN  @Path_in  TO  * 

1770  SUBEND 
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1000  !  PROGRAM:  PLOT_PUS 

1005  ! 

10)0  i  -  Routine  plots  drifter’s  position  and  u-component  speed. 

1015  i 

,02®  i  -  INPUT  DATA  SECTION - 

I0Z5  i  Title®:  Title  of  this  plot 
1030  DATA  "Drifter  Position  &  U_speed" 

1035  i  - 

1040  i  X_label$:  X_axis  Label 
1045  DATA  "X-DIST.  <N)" 

1050  i  - 

1055  i  Y_1 abel $ :  Y_axis  label 
1060  DATA  "Y-DIST.  <M)" 

1065  i  - 

1070  i  X_lef t , X_r ight , X_t ick_spacing , X_maj or_count 
1075  DATA  0,  4000,  Z50,  Z 

1 080  i  - 

1085  i  Y_bottom,  Y_top , Y_t ick_spaci ng , Y_maj or_count 
1080  DATA  0,  4000,  Z50,  Z 

1095  i  - 

1100  !  Grid:  Need  grid  lines  < yes/noj  1/0) 

1 105  DATA  0 

1110  !  - 

1115  !  Dump_graphics :  Dump  graphics  to  printer  (yes/no;  1/0) 

11Z0  DATA  1 

11ZS  !  - 

1130  I  Xoff.Yoff?  offset  ir  the  x-  and  y-  directions  respectively. 

1135  DATA  Z 34500, 349000 

1140  !  - - 

1145  !  Date_$:  plotting  date 

1150  DATA  1  APR  1985 

1155  i  - - - 

1160  1  Time_$:  plotting  tine 

1 165  DATA  14:30:00 

1170  i  - 

1 175  )  Ndrifter:  number  of  drifters  to  be  plotted,  and  titles  of  the  drifters 

1  180  DATA  6,  5,  6,  7,  8,  9,  10 

1185  i  - 

1190  !  Input_file$:  input  file  name 

1195  DATA  "DRIFTER_5= .700, 1” 

1196  DATA  "DRIFTER_6: ,700, 1" 

1Z00  DATA  "0RIFTERJ7: ,700, 1" 

1 Z05  DATA  "0RIFTER_8: ,700, I" 

1  Z 1 0  DATA  "0RIFTER_9:  ,700, 1" 

1 Z 1 1  DATA  "ORIFTER_I0: ,700, t” 

1  Z  1  5  ! 

IZZ0  i  - END  OF  INPUT  DATA - 

1 ZZ5  ! 

1  Z 30  OPTION  BASE  1 
1  Z 35  DIM  Tit  1 e$1 301 , 1 d*I Z ] 
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1240  DIM  Date Jt  1  1  ]  ,Time_$I8] 

1245  READ  Titles 
1Z50  READ  X_label$ 

1255  READ  Y_label$ 

12G0  REAO  X_1  ef  t  ,  X_r  ight .  X_t  ick_spacing  ,  X_maj  or_count 

1265  READ  Y_bottom , Y_top , Y_t ick_spacing , Y_maj or_count 

1270  READ  Grid 

1275  READ  Dump_graph i cs 

1280  READ  Xof f , Yof f 

1285  READ  Date_$ 

1290  READ  T ine_$ 

1295  Time=DATE<  Date_$  >  +  TIME< Time_$> 

1300  Major_tick_size=3.0 
1305  C$=CHR$( 255)&"K" 

1310  OUTPUT  Z  USING  I  Clear  screen  for  graph 

1315  GINIT  !  Initialize  various  graphics  parameters . 

1320  !  PLOTTER  IS  3 INTERNAL"  !  Use  the  internal  screen 

1321  PLOTTER  IS  705,"HPGL"  !  Use  the  external  HPGL  pen  plotter  at  HP-IP 

Address  705 

1325  GRAPHICS  ON  !  Turn  on  the  graphics  screen 

1330  LORG  6  !  Reference  point:  center  of  top  of  label 

1335  !  Determine  hou  many  GDUs  uide  and  high  the  screen  is 

1 340  Gdu<  X_gdu_max , Y_gdu_max ) 

1345  FOR  I=-.1  TO  .1  STEP  .1  !  Offset  of  X  from  starting  point 

1350  MOVE  X_gdu_max/2 . 3+1 , Y_gdu_max  !  Move  to  about  middle  of  top  of  screen 

1355  LABEL  USING  "#,K" ; Title*  !  Write  title  of  plot 

1360  NEXT  I  I  Next  position  for  title 

1365  DEG  !  Angular  mode  is  degrees  (used  in  LDIR) 

1370  Label! 3. 8, . 6,90,6, 1 ,0. ,Y_gdu_max/Z ,Y_label$)  !  Write  Y_axis  label 
1375  Label ( 3. 8 , . 6 ,0,4 , 1 , X_gdu_max/Z , .01 *Y_gdu_max , X_label$>  !  Write  X_axis  labe 
1 

1 380  Label (3. 8, .6, 0,4,1 ,. 85*X_gdu_max , . 8*Y_gdu_max , Date_$ ) 

1 385  Label (3. 8, .6, 0,4, I,. 85*X_gdu_max , . 75* Y_gdu_max , Time_$ ) 

1390  Label (3. 8, .6, 0,4,1 ,. 87*  X_gdu_max , . 7*  Y_gdu_max , " X_of  f : " &VAL$<  Xof  f  > ) 

1  395  Label  (3. 8, .6, 0,4,1,. 87*X_gdu_max , . 65*Y_gdu_max , M Y_of f : "&VAL$<  Yof f ) ) 

1400  I  Define  subset  of  screen  area 

1405  VIEWPORT  . 1 Z*X_gdu_max , .75*X_gdu_max , . 1 *Y_gdu_max , . 93*Y_gdu_max 
1410  l  Anisotropic  scaling:  left/right/bottom/top 
1415  WINDOW  X_lef t ,X_right , Y_bottom , Y_top 
14Z0  !  Draw  a  box 

14Z5  AXES  X_tick_spacing,Y_tick_spacing,X_left ,  Y_bottom ,  X_maj  or_r.ount  , Y_major_c 
ount , Ma’ ^r_t ick_s i ze 

1 430  AXjiS  X_t  ick_s pacing ,  Y_t  ick_spacing ,  X_right ,  Y_top ,  X_maj  or_count ,  Y_maj  or_cou 
nt , Maj  or_t ick_size 

1435  IF  Grid  THEN  GRID  X_t ick_spacing , Y_t ick_spacing , X_1 ef t , Y_bottom , X_maj or_co 
unt , Y_maj  or_count 

1440  CLIP  OFF  !  So  labels  can  be  outside  VIEWPORT  limits 

1445  CSIZE  3.5, .6  !  Smaller  chars  for  axis  labelling 

1450  LORG  6  !  Ref.  pt :  Top  center  !\ 

1455  WINDOW  X_lef t , X_right , . 1 *Y_gdu_max , . 9*Y_gdu_max 
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1460  X_step=X  t ick  _spacing*X  jiaj or_ccunt 

1465  FOR  I =X_ l ef  t  TO  X_nght  STEP  X_step'  Every  X_STEP  units 
1470  MOVE  I  -  X _1 ef t + I  , , 09  * Y_gdu_max 1  ft  smidgeon  belou  X-axis 
1475  LftBEL  USING  "#,K" ; I  I 

1480  NEXT  I  i 

1485  WINDOW  X_1 ef  t , X_r lght , Y_bottom ,  Y  top 

1490  LORG  8  I  Ref.  pt :  Right  center 

1495  Y_step=Y_t ick  spacing* Y_maj or_count 
1500  FOR  I=Y_bottom  TO  Y_top  STEP  Y_step  ! 

1505  MOVE  .8,1  1  Smidgeon  left  of  Y-axis 

1510  LftBEL  USING  "K" ; VAL$( I  )&"  "  I 
1515  NEXT  I  I 

1 5Z0  PENUP 

15Z5  1  ftnisotropic  scaling:  I ef t /r ight/bott om/ t op 
1530  WINDOW  X_1 ef  t , X_r ight , Y_bottom , Y_  top 
1535  ! 

1540  !  -  get  data  to  be  plotted 

1545  ! 

1550  REftD  Ndrifter 

1551  ftLLOCflTE  INTEGER  Dr  if ters( Ndrifter > 

1  55Z  REftD  Drifters! *> 

1555  FOR  Icurve=1  TO  Ndrifter 
1560  REftD  Input_f ile$ 

1 565  Get_data(  Input_f ile$ , Id$ , X_pos , Y_pos , Speed ,Xoff,Yoff,Time, Ndata) 

1570  IF  Ndata  THEN 

1571  Id$=VAL$< Drifters!  Icurve ) ) 

1575  Label! 3.8, .6,0,4, 1 , X_pos , Y_pos+ 1 50 , Id$ > 

1580  Label!  3.8, .6,0,4, 1 , X_pos , Y_pos- 1 50.  ,VAL$< INT! Speed* 100)/ 100) ) 

1585  Label! 3, .6,0,4, 1 , X_pos , Y_pos , "  +  “ ) 

1590  PENUP 

1595  END  IF 
1600  NEXT  Icurve 
1 605  BEEP 

1610  ON  KBD  GOTO  Exit 

1615  DISP  "Enter  'Space  bar’  to  go  on" 

I6Z0  Idle:G0T0  Idle  1  View  the  plot  as  long  as  you  want 

1625  Ex  it : DI SP 

1630  IF  Dump_graphics  THEN 

1635  !  DUMP  GRAPHICS  CRT  TO  #701 

1640  END  IF 

1645  GRAPHICS  OFF 

1650  OUTPUT  Z  USING  "#,K";C$ 

1655  END 
1660  ' 

1665  Gdu:SUB  Gdu!  X_gdu_max ,Y_gdu_max .OPTIONAL  Gdu_xmid ,Gdu_ymid ) 

1670  !  This  returns  Xright,  Yhigh  and  their  respective  midpoints  in  GDUs . 

1675  1  Note  that  if  Gud_xmid  is  defined,  Gdu_ymid  must  be  also. 

1680  COM  /G_units/  Gdu_xmax ,Gdu_ymax , Udu_xmi n , Udu_xmax , Udu_ymi n , Udu_ymax , Show 

1685  IF  Gdu_xmax=0  THEN 

1690  Gdu_xmax=100*MAX< 1 .RATIO) 


!  \ 

!  >  Label  X-axis 

!  / 

1  / 

!  \ 

!  \ 

1  >  Label  Y-axis 

!  / 

1  / 
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1 695  Gdu_ymax=100*MAX( 1 , I ./RATIO) 

1700  END  IF 
1705  X_gdu_max =Gdu  _.xnax 
1710  Y_gdu_max=Gdu_ymax 
1715  IF  NPAR>2  THEN 
1720  Gud_xmid=Gdu_xmax* . 5 

1725  Gud_ymid=Gdu_ymax * . 5 

1730  END  IF 
1735  SUBEND 

1740  Label : SUB  Label ( Os ize , Asp_rat io , Ld ir , Lorg , Pen , X , Y , Tex t$ ) 

1745  1  This  defines  several  systems  variables  (in  CSIZE,  LDIR,  etc.),  and 

1750  1  labels  the  test  (if  any)  accordingly. 

1755  DEG 

17G0  CSIZE  Csize,Asp_ratio 
1 7G5  LDIR  Ldir 

1770  LORG  Long 

1775  PEN  Pen 
1780  MOVE  X ,  Y 

1785  IF  Text$<>""  THEN  LABEL  USING  ;  Texts 

1730  PENUP 
1795  SUBEND 
1800  ! 

1 805  Get_data: SUB  Get_data( File_in$ , Drno$ , X_pos , Y_pos , Speed , Xof f , Yof f , Time .Ndata 
) 

1810  ! 

1815  OPTION  BASE  1 
1820  DIM  Id*[ 1 1 ,Oesc$[ 30] 

1 8Z5  INTEGER  Npos 

1830  ASSIGN  @Path_in  TO  File_in$ 

1835  ENTER  @Path_in ; Drno$ , IdS , Oesc$ , Npos 
1840  Ndata=0 

1845  ON  END  @Path_in  GOTO  Close_file 

1850  ENTER  @Path_in; X_pos , Y_pos , T , Speed , Dummy  I  Change  Dummy  and  Speed 

for  U&V 

1855  IF  TOTime  THEN  1845 
1 8G0  X_pos=X_pos-Xof f 
1865  Y_pos=Y_pos-Yoff 
1870  Ndata= 1 

1875  Close_file: ASSIGN  @Path_in  TO  * 

1880  SUBEND 
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1000  i  PROGRAM:  PLOT_PVS 
1005  i 

1010  i  -  Routine  plots  drifter’s  position  and  u-component  speed. 

1015  i 

10Z0  i  -  INPUT  DATA  SECTION - 

I0Z5  i  Title®:  Title  of  this  plot 
1030  DATA  "Drifter  Position  li  V__speed“ 

1035  i  - 

1040  I  X_label$:  X_axis  Label 
1045  DATA  "X-DI5T.  (  N ) " 

1050  i  - 

1055  i  Ylabel®:  Y_axis  label 
1060  DATA  "Y-OIST.  (M)M 

1065  i  - 

1070  !  X_1 ef t , X_r lght , X_t ick_spaci ng , X_maj or_count  • 

1075  DATA  0.  4000,  Z50,  Z 

1080  i  - 

1085  !  Y_bottom,  Y_t op , Y_t ick_spac i ng , Y_maj or_count 

1090  DATA  0,  4000,  Z50,  Z 

1095  i  - 

1100  1  Grid:  Need  grid  lines  (yes/no;  1/0) 

1 105  DATA  0 

1110  i  - - 

1115  I  Dump_graphics :  Dump  graphics  to  printer  (yes/no;  1/0) 

1 1Z0  DATA  1 

11ZS  !  - - - 

1130  i  Xoff.Yoff?  offset  in  the  x-  and  y-  directions  respectively. 
1135  DATA  Z 34500, 349000 

1140  !  - 

1145  !  Date_$:  plotting  date 

1150  DATA  1  APR  1985 

n55  i  - 

1160  1  Time_$:  plotting  time 

1 165  DATA  14: 30:00 

1170  !  - - - 

1175  1  Ndrifter:  number  of  drifter  to  be  plotted 

1 1 80  DATA  4 

,185  i  - - 

1190  1  Input_file$:  input  file  name 

1195  DATA  "DRIFTER_5: ,700, 1” 

IZ00  DATA  ''0RIFTER_7:  ,700,  1” 

1  Z05  DATA  ’‘DRIFTER_8:  ,700,  1" 

IZ10  DATA  "DRIFTER_9: ,700, I" 

1Z15  I 

1ZZ0  i  - END  OF  INPUT  DATA - 

I ZZ5  i 

1 Z 30  OPTION  BASE  1 
I Z 35  DIM  Ti 1 1 e$[ 30] ,Id$(Z] 

1Z40  DIM  Date_$t  1  1  1  , Time_$t 8] 

1 Z45  READ  Title® 
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1250  READ  X_1 abel $ 

1255  READ  Y_label$ 

1260  READ  X_1  ef  t  ,  X_r  ight ,  X_t  ick__spacing ,  X_maj  or__count 

1265  READ  Y_bottom  ,  Y_top  ,  Y_t  ick_spacing  ,  Y_.maj  or_count 

1270  READ  Grid 

1275  READ  Dump_graphi cs 

1280  READ  Xoff .Yoff 

1285  READ  Date_$ 

1290  READ  TimejS 

1295  T ime=DATE<  Date_$ )+TIME( Time_$) 

1300  Maj or_t ick_s ize=3 . 0 
1305  C$=CHR$(  255  )&"K" 

1310  OUTPUT  2  USING  "#,K";C$  i  Clear  screen  for  graph 

1315  GINIT  !  Initialize  various  graphics  parameters. 

1320  PLOTTER  IS  3 ," INTERNAL"  !  Use  the  internal  screen 

1325  GRAPHICS  ON  !  Turn  on  the  graphics  screen 

1330  LORG  6  !  Reference  point:  center  of  top  of  label 

1335  !  Determine  how  many  GOUs  wide  and  high  the  screen  is 

I  340  Gdu(  X_gdu_max , Y_gdu_max ) 

1345  FOR  I=-.2  TO  .2  STEP  .1  !  Offset  of  X  from  starting  point 

1350  MOVE  X_gdu_max/2 . 3+1 , Y_gdu_max  !  Move  to  about  middle  of  top  of  screen 

1355  LABEL  USING  "# , K" ; Ti t le$  !  Write  title  of  plot 

1360  NEXT  I  !  Next  position  r title 

1365  DEG  !  Angular  mode  .  iegrees  (used  in  LDIR) 

1370  Label ( 3. 8, . 6, 90, 6, 1 , 0. , Y_gdu_max/2 , Y_label $)  ;  Write  Y_axis  label 

1375  Label ( 3. 8 , . 6 ,0,4 , 1 ,X_gdu_max/2 , .01 *Y_gdu_max , X_label$ )  !  Write  X_axis  labe 
1 

1 380  Label <  3.8, .6.0,4, 1 ,. 85*X_gdu_max , . 8*Y_gdu_max ,Date_$) 

1 385  Label (3. 8, .6, 0,4,1,. 85*X_gdu_max , . 75* Y_gdu_max , Time_$ ) 

1 390  Label ( 3.8, .6,0,4, 1 , . 87*X_gdu_max , .7*Y_gdu_max , "X_of f : "&VAL3K  Xoff ) > 

1  395  Label(  3.8,  .6,0,4, 1  ,  .  87*X__gdu_max  , .  65*Y_gdu_max  ,  "  Y_of  f  :  "&VALS(  Yoff  >  ) 

1400  !  Define  subset  of  screen  area 

1 405  VIEWPORT  . 1 2*X_gdu_max , . 75*X_gdu_max , . 1  * Y_gdu_max , . 9  3* Y_gdu_max 
1410  1  Anisotropic  scaling:  lef t/right/bottom/top 
1415  WINDOW  X_lef t ,X_right , Y_bottom , Y_top 
1420  !  Draw  a  box 

1425  AXES  X_t ick_spacing , Y_t ick_spacing,X_lef t , Y_bottom , X_maj or_count ,Y_major_c 
ount , Maj or_t ick_5ize 

1 430  AXES  X_t ick_spacing , Y_t ick_spacing , X_r ight , Y_top , X_maj  or_count , Y_maj  or_cou 
nt , Maj  or_t ick_size 

1435  IF  Grid  THEN  GRID  X_t ick_spacing , Y_t ick_spacing , X_1 ef t , Y_bottom , X_maj or_co 
unt , Y_maj  or_count 

1440  CLIP  OFF  !  So  labels  can  be  outside  VIEWPORT  limits 

1445  CSIZE  3. 5,. 6  !  Smaller  chars  for  axis  labelling 

1450  LORG  6  !  Ref.  pt:  Top  center  !\ 

1455  WINDOW  X_lef t , X_r ight , . 1  * Y_gdu_max , . 9* Y_gdu_max 
1 460  X_step=X_t ick_spacing*X_maj  or_count 

1465  FOR  I=X_lef t  TO  X_right  STEP  X_step!  Every  X_STEP  units  1  \ 

1470  MOVE  I -X_lef t+ 1 , . 09* Y_gdu_max !  A  smidgeon  below  X-axis  I  >  Label  X-axis 

1475  LABEL  USING  "#,K";I  !  !  / 
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1480  NEXT  I  i  ,'/ 

1485  UINDOU  X_1 ef t , X_r lght . Y_bot t on , Y_top 

1490  LORG  8  !  Ref.  pt:  Right  center  !\ 

1495  V_step=Y_t iek_spaci ng* Y_maj or_count 

1500  FOR  I = Y_bottom  TO  Y_top  STEP  Y_step  '  !  \ 

1505  MOVE  .8,1  !  Smidgeon  left  of  Y-axis  1  >  Label  Y-axis 

1510  LABEL  USING  "K" ; VAL$( I  )&"  '  !  / 

1515  NEXT  I  i  I  / 

1 5Z0  PENUP 

1525  1  Anisotropic  scaling:  1 ef t / r ight/bot tom/ top 
1530  UINDOU  X_1 ef t , X_r lght , Y_bot t on , Y_t op 
1535  i 

1540  1  -  get  data  to  be  plotted 

1545  ' 

1550  READ  Ndnfter 

1555  FOR  I  cur ve=  1  TO  Ndnfter 

15G0  READ  I nput_f i 1 e$ 

1  5G5  Get_data(  Input_f ile$ , Id$ , X_pos , Y_pos , Speed, Xoff ,Yoff .Time, Ndata ) 

1570  IF  Ndata  THEN 

1575  Label <  3.8, .6,0,4, 1 , X_pos , Y_pos+ 1 50 , I d$ > 

1580  Label ( 3. 8,. G, 0,4,1 , X_pos , Y_pos- 1 50. ,VAL$< INT( Speed* 100)/ 100) ) 

1585  Label ( 3 , ,6,0,4, 1 , X_pos , Y_pos , " +“ ) 

1590  PENUP 

1595  END  IF 
1 G00  NEXT  Icurve 
1605  BEEP 

1610  ON  KBD  GOTO  Exit 

1615  DISP  "Enter  ’Space  bar’  to  go  on" 

1620  Idle:GOTO  Idle  !  Vieu  the  plot  as  long  as  you  want 

1  6Z5  Ex  it : DI 5P 

1630  IF  Oump_graphics  THEN 

1635  DUMP  GRAPHICS  CRT  TO  #701 

1640  END  IF 

1645  GRAPHICS  OFF 

1650  OUTPUT  2  USING  "#,K" ;C$ 

1655  END 
1660  ! 

1665  Gdu:SUB  Gdu< X_gdu_max , Y_gdu_max , OPTI ONAL  Gdu_xmid,Gdu_ymid) 

1670  1  This  returns  Xright,  Yhigh  and  their  respective  midpoints  in  GDUs. 

1675  !  Note  that  if  Gud_xmid  is  defined,  Gdu_ymid  must  be  also. 

1680  COM  /G_units/  Gdu_xmax ,Gdu_ymax , Udu_xmin , Udu_xmax , Udu_ymin , Udu_ymax , Show 

1685  IF  Gdu_xmax=0  THEN 

1690  Gdu_xmax= 1 00* MAX< 1 .RATIO) 

1695  Gdu_ymax= !00*MAX(  1,1. /RATIO ) 

1700  END  IF 

1705  X_gdu_max=Gdu_xmax 

1710  Y_gdu_max=Gdu_ymax 

1715  IF  NPAR>2  THEN 

1720  Gud_xmid=Gdu_xmax* . 5 

1725  Gud_ymid=Gdu_ymax» . 5 
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1730  END  IF 
1735  SUBEND 

1740  Label:SUB  Label  (  Cs lze , Asp_rat 10 ,Ldir , Lorg , Pen , X , V , Text® ) 

1745  I  This  defines  several  systems  variables  (in  CSIZE,  LDIR,  etc.),  and 
1750  1  labels  the  test  (if  any)  accordingly. 

1755  DEG 

I7S0  CSIZE  Cs iz e , Asp_rat l o 
1 7G5  LDIR  Ld i r 

1770  LORG  Lorg 

1775  PEN  Pen 
1780  MOVE  X  ,  Y 

1785  IF  Text$<>“"  THEN  LABEL  USING  ;  Text$ 

1790  PENUP 
1795  SUBEND 
1 800  i 

I  805  Get_dat a: SUB  Get_data(  F i 1 e_i n$ , Drno$ , X_pos , Y_pos , Speed ,Xoff,Yoff,Ti me , Ndata 
) 

1810  i 

1815  OPTION  BASE  1 

1 8Z0  DIM  Id$I 1 1 , Desc$I 301 

1825  INTEGER  Npos 

1830  ASSIGN  @Path_in  TO  File_in$ 

1835  ENTER  @Path_in; Drno$ , I d$ , Desc$ , Npos 
1840  Ndata=0 

1845  ON  END  @Path_in  GOTO  Closest ile 

1850  ENTER  @Path_in ; X_pos , Y_pos , T , Dummy , Speed 

1855  IF  TOTime  THEN  1845 

I860  X_pos=X_pos-Xof f 

1865  Y_pos=Y_pos- Yof f 

1870  Ndata= 1 

1875  Close_f  ile:  ASSIGN  @Path__in  TO  * 

1880  SUBEND 
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1000  i  Program  PREPARE 

1005  DIM  Id$[ 1 1 ,Desc$[ 30] ,Drid$I Z ] 

1010  INTEGER  Npos . I .Flag, Check ,N 
1015  REAL  X  ,  Y  ,  T  ,  T 1 ,U,V 
10Z0  DIM  Length!  14) 

1025  DATA  28,45,64,38,51,99,100,78,139,159,154,186,165,47 
1030  ASSIGN  @Path_in  TO  “LDI01 APR: ,700, 1" 

1035  FOR  I dr= 1  TO  14 
1040  Dnd$=VAL$<  Idr) 

1045  READ  Length! Idr  ) 

1050  CREATE  BOAT  "DRIFTER_',&7AL$(  Idr  )&“  :  ,700,  1  "  .Length!  Idr  )  ,80 
1055  ASSIGN  @Pa th_out  TO  "DRI FTER_"&VAL$( Idr )&": ,700, 1 " 

1060  ENTER  @Path_in; Id$ ,Desc$ ,Npos 

1065  PRINT  Drid$,Id$,Desc$,Npos 

1070  OUTPUT  @Path_out j IdridS , Id$ , Desc$ , Npos 

1075  ENTER  @Path_in; FI ag , X 1 , Y 1 , T 1 

1080  ENTER  @Path_in; Flag , X , Y , T 

1085  U=< X-X1 )/120,*100. 

1090  V=< Y-Y1  >/1Z0.*100. 

1095  X1=X 

1100  Y 1 =Y 

1105  T 1 =T 

1 1 '0  PRINT  X , T , U , V 

1115  OUTPUT  @Path_out ; X , Y , T ,U , V 

1 120  FOR  1=3  TO  Npos 

1125  ENTER  @Path_in; Flag, X.Y.T 

1130  U=!X-X1 ) / <  T-T 1  >*100. 

1135  V=< Y-Y1 )/<  T-T1  )*100. 

1140  X 1 =X 

1145  Y 1  =  Y 

1150  T 1  =T 

1  1 55  OUTPUT  @Path_out ; X , Y , T , U , V 
1  160  PRINT  I  ,X,T,U 

1 165  NEXT  I 

1 170  Next_idr : NEXT  Idr 
1175  END 
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1000  I  PROGRAM:  PRINTJDA 

1001  I  -  routne  prints  the  Objective  analysis  results. 

1005  OPTION  BASE  I 

1010  DIM  Ti  1 1  e  1  $[  801  , Tit le2$l  801  .  A__$I  25B1  ,  App_$L  2561 

1015  DIM  Phi( Z1 ,2 1 ) ,A( 15. 15) 

1016  Clear$=CHR$( 255)&"K" 

1017  OUTPUT  2  USING  "#,K" ; Cl  ear $  !  clear  the  screen. 

1020  i  - 

1025  !  assign  the  objective  analysis  results  file  to  Oa_file$ 

1030  INPUT  "Enter  the  filename  of  objective  analysis  file  to  be  pr i nted" , Oa_f i 1 

e$ 

1035  Oa_f ils$=UPC$< Oa_f ile$) 

1 040  !  - 

1045  I  read  in  the  x-  and  y-gnd  points  of  the  0.  A.  field 

1050  INPUT  "Enter  the  x-  and  y-grid  points  of  the  computing  f ield( NPX , NPY ) " , Npx 

.  Npy 

1055  IF  Npx <  =  Z  1  AND  Npy<'  =  21  THEN  1085 
1060  BEEP 

1065  PRINT  TABXY( 10, 10) *******  Warning  *******" 

1070  PRINT  TABXY( 1 0, 1 1 ) , "The  sizes  of  array  Phi(i.j)  defined  in  this  version" 
1075  PRINT  TABXY( 1 0 , 1 2  )  , "are  incorrect,  make  correction  in  lines  1015  and  1055. 

1080  PRINT  TABXY( 10, 1 3) ."Program  execution  is  halted." 

1085  STOP 

1080  !  set  the  scaling  factor  for  o.a.  output 

1035  INPUT  “Enter  the  scaling  factor  for  entire  field  data" , Seal ing_f actor 
1100  1  need  a  hard  copy7  (yes/no;  1/0) 

1105  INPUT  "Need  a  hard  copy?  (yes/no?  1/0>",Copy 
1110  IF  Copy  THEN  PRINTER  IS  701 

1115  !  - 

1120  ASSIGN  @Path_in  TO  Oa_file$ 

1125  ON  END  @Path_in  GOTO  Close_file 
1130  ENTER  @Path_in; Title! $ 

1135  ENTER  @Path_in; TitleZS 
1 140  FOR  I =Npy  TO  1  STEP  -1 
1 145  FOR  J=1  TO  Npx 
1150  ENTER  @Path_in;Phi< I , J) 

1 155  NEXT  J 
1 160  NEXT  I 

1165  CALL  Intrp(Phi( * ) , 15, A<  * ) ) 

1170  OUTPUT  2  USING  ”#,K";C$ 

1175  PRINT  TAB< 30) ,Title1$ 

1180  PRINT  TAB(  30) ,T’tle2$ 

1185  PRINT  USING  "K , SO. OE" ; "Seal e : " , Seal i ng_f act or 

1180  PRINT  "  1  2  3  4  5  6  7  8  3  10  1  1  12  1 

3  14  15" 

1  135  PRINT  "  ************ 

*  #  #  ” 

1200  FOR  1=1  TO  15 
1205  A  $=VAL$(  16-1  )&"*" 
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1210  IF  LEN(  ft  $  >  3  I  HEN  "&AJE 

1215  FOR  J= I  TO  15 

1  220  Rpp VHL$(  INK  rU  I  , J ) / Seal l nq_f act  or  )  ) 

1225  SELECT  LEN<  App._$  > 

1230  CASE  0 

1235  CASE  1 

1240  App_$="  "iiHpp  $ 

1245  CASE  2 

1250  App_$="  "&App_$ 

1255  CASE  3 

1260  App_$="  "iApp_$ 

1265  CASE  4 

1270  App_$="  " &App  _$ 

1275  END  SELECT 

1280  A_$=A__$&App_$ 

1285  NEXT  J 
1280  PRINT  A_$ 

1285  NEXT  I 

1300  DISP  "Enter  ’ CONTONUE’  to  go  on" 

1 305  PAUSE 
1310  DISP 
1315  GOTO  1 125 

1320  C 1 ose_f l 1 e : ASSIGN  @Path_in  TO  * 

1325  PRINTER  IS  CRT 
1330  END 
1  335  l 
1 340  i 
1  345  ' 

1350  I ntrp : SUB  I n lrp< B< * ) , M , A< * ) ) 

1355  OPTION  BASE  1 

1360  INTEGER  Lli.Llipl ,Llj ,Lljp1 

1365  1  Routmne  interp'’ates  array  B(M,M>  onto  array  A(  15,15). 

1370  FOR  1-1  TO  15 

1375  LI i  =  INT( ( M- I )*< I- 1 )/ 14. +.38398  ) 

1380  LI  ip  I =L 1 1  +  1 

1 385  IF  LI i=0  THEN  L 1 x  = 1 

1390  DiL=( 1-1  ) *  <  M- 1  )/14.-(Lli-1  ) 

1385  Dir= 1 .-Dil 

1400  FOR  J=1  TO  15 

1405  LI j  =INT<  <  M- 1 >»< J-1 )/ 14.  +  .99333  ) 

1410  LIjp1=Llj+1 

1415  IF  L 1 j  =0  THEN  LI  j  =  1 

1420  Dj I =< J-!  >  *  (  M-  1  )/14.-(Llj-1 ) 

1425  Dj  u= 1 . - D j 1 

14  30  A(  I  ,1)  =  '.  Dxr»6(Lli,Llj  )+Dil*B(Llip1  ,Llj  ))*Dju  +  (Dir»B(Lli,Lljp1  )+Dil*B(Ll 

i  p  1  ,  L  1  j  p  I  )  >  *  Dj  1 

1435  NEXT  J 

1440  NEXT  I 

1445  SUBEND 
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1000  i  PROGRAM  GET_IP 

1005  i  ROUTINE  CALCULATES  INTERPOLATION  POSITION 
1010  i 

1015  OPTION  BASE  1 
I0Z0  RAD 
1025  Xof f  =  Z  34500 
1030  Yof  f  =  349000 
1035  'Input  dataport 
1070  I 

1075  '  -  create  a  BDAT  file  for  output 

1076  i  out put_f i 1 ename , f i 1 e_l ength , record_s i 2 e 

12-80  READ  Filename®  ,File_length  ,Record_size 

1081  DATA  "IP_P0S: ,700, I" ,225,80 

1082  ON  ERROR  GOTO  1084 

1083  CREATE  BDAT  F i 1 ename® , F i 1 e_l ength , Record_s l z e 

1084  OFF  ERROR 

1085  ASSIGN  @Path_out  TO  Filename® 

1090  i 

1095  READ  Npx.Npy  !  read  in  the  number  of  grid  points  in  x-  and  y-directions 

1096  DATA  15,15 

1  100  FOR  J=1  TO  Npy 
1105  Y=(  J- 1 )*Z50.  +  Yoff 

1  1 10  FOR  1  =  1  TO  Npx 
1115  X=< 1-1 )*Z50.+Xoff 

1  1  20  OUTPUT  @Path_out ; X , Y 
1  125  PRINT  X , Y 

1  1  30  NEXT  I 
1 1 35  NEXT  J 

1140  ASSIGN  @Path_out  TO 
1145  END 
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!  000  i  PROGRAM  SOA 
1005  i 

1010  1  »*****+*»*»»***************»***#»»**»*»*#***»****»**»*»if*»*****##*»»»* 


1015  i  *  * 
10Z0  1  *  Scalar  Space- Time  Objective  Analysis  Package  * 
1 0Z5  i  *  * 
1030  I  *  Language:  BASIC  3.0  » 
1035  1  »  System:  Hewlett  Packard  98IS  * 
1040  i  *  Version:  1.00  * 
1045  i  *  Date:  July  1986  * 
1 050  l  *  * 
1055  1  *  Revised  by  Dr.  L.  Charles  Sun  * 
1060  l  *  Hawaii  Institute  of  Geophysics  * 
1065  l  *  University  of  Hawaii  at  Manoa  * 
1070  !  *  Honolulu,  Hawaii  968Z2  * 
1 075  '  *  * 
1080  !  *  Developed  at  U.  S.  Coast  Guard  Research  and  Development  Center  * 
1085  I  *  Groton,  Connecticut  06340  * 
1090  I  *  * 


1095  !  ********************************************************************** 

1100  i 
1105  ! 

1110  OPTION  BASE  I 
1115  RAD 

1120  DIM  MessageSC 256] 

1125  DIM  Oate0$! 15] ,Time0$t 15] 

1130  DIM  Xdata! 2000) ,Ydata(  2000) ,Tdata< 2000) ,Phidata< 2000) 

1135  DIM  Xrd!  1089), Yrd!  1089), Trd( 1089 ) ,Phird( 1089) 

1140  DIM  Xopd! 1089 > ,Yopd( 1089 ) ,Topd< 1089 > .Phiopd! 1089 > 

1145  DIM  Xip<  1089) ,Yip( 1089) , Tip!  1089) 

1150  DIM  Theta! 1089 ) ,Es( 1089 ) ,Erms( 1089) 

1155  DIM  Cor!  1089 ), Index!  1089 ) 

1160  COM  /Fold/  Xfold,Yfold,Tfoid 
1 165  COM  /Ef ield/  Evar 
1 170  COM  /Cblock/  C<40> 

1175  i 

I  I 80  ClearS^CHRS!  255  )&"K" 

1185  SdateS= DATES! TIME DATE ) 

1190  St ime$= TIMES!  TIMEDATE ) 

1195  I 

1200  OUTPUT  KB0;Clear$;  1  clear  the  screen. 

1205  ' 

1210  PRINT  TABXY!  1 5 , 10) , "Do  you  need  documentations?" 

1215  INPUT  "Enter  Y/N  for  yes/no", AnsS 

1220  OUTPUT  KBD;Clear$;  !  clear  the  screen. 

1225  IF  UPCS! AnsSI 1 ; 1 ] )<>"Y"  THEN  GOTO  1450 
1 Z  30  Ty=  1  0 

1235  PRINT  TABXY!  1 0,7 ), "Which  section  do  you  want  to  look  at  ?" 

1240  PRINT  TAB!  Ty) ."Section  1:  INTRODUCTION" 

1245  PRINT  TAB!  Ty ) , "  2:  BASIC  THEORY" 
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1250  PRINT  TflB(Ty),"  3:  DETERMINATION  OF  CORRELATON  FUNCTION- 

1255  PRINT  TAB(Ty>,"  4:  ELEMINATION  OF  DISTANT  DATA" 

1260  PRINT  TAB(  ty )  , "  5:  PROGRAMS  DESCRIPTIONS" 

1 2G5  PRINT  TAB(Iy),"  6:  USER  INSTRUCTIONS  AND  EXAMPLES" 

1270  PRINT  TAB(Ty),"  7:  REFERENCES" 

1275  PRINT  TAB<  Ty ) , "  8:  DEFINITIONS  OF  VARIABLES" 

1280  PRINT 

1285  PRINT  IAB( Ty ), "Enter  E  to  exit  help  nanus ,  and  start  computation." 

1290  INPUT  " Sel ect ion? ",  Code® 

1295  OUTPUT  KBD; C 1  ear® i 

1300  IF  Code$="E"  OR  Code$="e"  THEN  GOTO  1450 
1305  F i le$="Doc_"&Code$ 

1310  ASSIGN  @P,ead_message  TO  File$ 

1315  Read_message:  ! 

1 320  Count  =0 

1325  ON  END  @Read_message  GOTO  End_mess 
1330  ENTER  @Read_message ; Message® 

1335  PRINT  TA8<7> , Message®! 8 ; 801 
1340  Count=Count+ I 
1345  IF  Count= 1 8  THEN 
1350  INPUT  "More  message  ?  ",Ans$ 

1355  IF  Ans$="y"  OR  Ans$="Y"  THEN 
1360  OUTPUT  KBD; Clear®; 

1365  GOTO  Read_message 

1370  ELSE 

1375  OUTPUT  KBD; Clear®; 

1380  GOTO  1235 

1385  END  IF 
1390  END  IF 
1395  GOTO  1325 
1400  End_mess:  1 

1405  INPUT  "End  of  this  section.  Need  another  section?  ",Ans$ 

1410  IF  Ans$=" y"  OR  Ans$="Y“  THEN 
1415  OUTPUT  KBD; Clear®; 

1420  GOTO  1235 
1425  ELSE 

1430  OUTPUT  KBD; Cl  ear®; 

1435  GOTO  1450 
1440  END  IF 
1445  I 

1450  I nput _paramet er  :  ! 

1455  OUTPUT  KBD; Clear®; 

1460  PRINT  "»**♦***»*»»»************»*************«************♦*•*♦*********»* 
*  #  *  " 

1465  PRINT  "* 

»  " 

1470  PRINT  "*  Scalar  Space-Time  Objective  Anal,  ackage 

*  " 

1475  PRINT  "* 

#  " 
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Language:  BASIC  3.01 


1480  PRINT  "* 

»  " 

1485  PRINT  "*  System:  Hewlett  Packard  9000-Z000 

*  '* 

1490  PRINT  "»  Version  1.00 

#  " 

1495  PRINT  "»  July  198G 

#  " 

1500  PRINT  "* 

*  “ 

1505  PRINT  ”*  Coast  Guard  Research  and  Development  Center 

*  " 

1510  PRINT  "*  Groton,  Connecticut  0G340 

*  '* 

1515  PRINT  "# 

*  " 

1 5Z0  PRINT  "***#********»**#***»»*****»******************#****************»***» 
*  *  #  " 

1 5Z5  PRINT 
1530  Tab= 1 0 

1535  PRINT  "Before  you  go  on,  please  check  the  followings:" 

1540  PRINT  "  1 )  Is  the  disc  containing  the  observed  data  and" 

1545  PRINT  "  the  inter/extrapolation  positions  in  DRIVE#  1  ?" 

1550  PRINT  "  Z)  Are  the  output  files  existed  ?,  if  not,  create  them  before" 
1555  PRINT  "  you  do  the  analysis." 

1 5G0  PRINT  ”  3)  The  file  length  for  each  output  should  be  Z  records  longer  tha 

n" 

1555  PRINT  "  the  total  inter/extrapolation  position  points." 

1570  PRINT  "  4)  The  current  version  allows  Z000  input  data  points  " 

1575  PRINT  "  and  1089  inter/extrapolation  points." 

1580  PRINT  ”  5)  This  version  gets  the  correlation  function  from  the" 

1585  PRINT  "  fitted  formula." 

1590  (PRINT 

1595  i PRINT  TAB( Tab) ."Send  your  inquires  to” 

1600  'PRINT  TAB<  Tab) , “  Dr.  L.  Charles  Sun" 

1505  'PRINT  TAB(Tab),"  Deptartment  of  Oceanography" 

1610  'PRINT  TAB( Tab ) , "  University  of  Hawaii" 

1515  'PRINT  TAB( Tab ) , "  Honolulu,  HI  9S8ZZ" 

1 6Z0  'PRINT  TAB( Tab ) , "  <808)948-7533" 

I6Z5  INPUT  "Hit  I  RETURN]  or  [ENTER]  to  cont i nue" , Answ$ 

1630  Tx  =  1 0 
1635  Ty= 1 0 
1640  ' 

1645  Lp_flg=0 
1650  LOOP 

1655  OUTPUT  KBD, •Clear®; 

1660  !  - 

1665  !  Date0$ , Time0$  =  time  of  the  analysis  to  make. 

1570  PRINT  TABXY< Tx ,Ty) ."Enter  time  of  the  analysis  to  make.(DD  MMM  YY , HH: MM: S 
S>" 
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1675  PRINT  IAB( ly ) , " ( end osed  with  quotations.)" 

1680  DISP  "Date0$? ,Time0$?" i 

1685  IF  Lp_f  lg  THEN  OUTPUT  KBD; . ;Date0$; . , . ;Time0$; . : 

1690  INPUT  "" ,Date0®,Time0® 

1695  Dats0$  =  UPC$(  Oate0$) 

1700  OUTPUT  KBD; Clear®; 

1 705  i  - 

1710  I  Nobj  =  number  of  objective  analysis  to  make. 

1715  PRINT  Tf)BXY(  Tx  ,  Ty ),  "Enter  the  number  of  objective  analyses  to  make." 

1720  DISP  "Nobj  ; 

1725  IF  Lp_flg  THEN  OUTPUT  KBD; Nobj ; 

1730  INPUT  ""  , Nobj 
1735  OUTPUT  :<8D;Clear$; 

1740  i  -  - 

1745  i  Delta_t  =  time  interval  for  each  extrapolation  in  time. 

1750  PRINT  fRBXY( Tx , Ty ). "Enter  the  time  interval  for  extrapolation  in  time" 

1755  PRINT  TflB(Ty),"  or  zero  for  instantaneous  computation." 

1760  DISP  "Delta_t?" ; 

1765  IF  Lp_f lg  THEN  OUTPUT  KBD;Delta_t; 

1770  INPUT  "" ,Delta_t 
1775  OUTPUT  KBD; Clear®; 

1780  i  - 

1785  i  Xlimit  =  Max.  distance  radius  from  the  reference  point  of  the  domain. 
1790  PRINT  TRBXY< Tx , Ty ), “Enter  maximum  distance  radius  from  the  reference  poin 
t  of  the  domain." 

1795  PRINT  TRB(Ty),"Rll  data  within  this  range  are  retained." 

1800  DISP  "Xlimit?"; 

1805  IF  Lp_f lg  THEN  OUTPUT  KBD;Xlimit; 

1810  INPUT  "" .Xlimit 
1815  OUTPUT  KBD; Clear®; 

1820  i  - 

1825  1  Tlimit  =*  Max.  time  radius  before  and  after  the  time  of  the  analysis. 

1830  PRINT  TRBXYI Tx , Ty ), "Enter  the  max.  time  radius  before  and  after  the  time 
of  the  analysis." 

1835  PRINT  TAB(Ty),"All  data  within  this  range  are  retained." 

1840  DISP  "Tlimit?"; 

1845  IF  Lp__f lg  THEN  OUTPUT  KBOiTlimit; 

1850  INPUT  "".Tlimit 
1855  OUTPUT  KBD; Clear®; 

1 860  I  - 

1865  !  Max_space_l ag  =  maximum  space  lag. 

1870  i  Max_time_lag  =  maximum  time  lag. 

1875  PRINT  TRBXY( Tx , Ty ) , " Set  the  influence  domain  for  computation  points," 

1880  PRINT  TR8(Ty),“The  maximum  time  lags  >=  the  multiplication  of  time  mterv 
al " 

1885  PRINT  TRB(Ty),"  and  the  number  of  the  analysis  to  make." 

1890  PRINT  TAB<Ty),"The  maximum  space  lags  >=  the  sizes  of  the  computation  dom 
ain .  ” 

1395  PRINT 

1900  PRINT  TRB( Ty ), "Enter  the  maximum  space  and  time  lags." 


B-28 


1905  DISP  " Max _space_l ag? , Max_ t ime _i ag? " ; 

1910  IF  Lp_flg  THEN  OUTPUT  KBDi Max_space_ 1 ag , Max_t ime_l ag ; 

1915  INPUT  "  "  ,  Max  __space_l  ag  ,  Max_t  ime_l  ag 
t  9Z0  OUTPUT  KBDi Clear®; 

1 9Z5  i  - 

1930  1  Limit  =  Max.  number  of  influential  points. 

1935  PRINT  TA8XY(  Tx ,  Ty ),"  Enter  the  maximum  number  of  influential  points  (40  ma 
x  ) .  " 

1940  DISP  ,,Limit?,‘; 

1945  IF  Lp_flg  THEN  OUTPUT  KBDjLimit; 

1950  INPUT  "" .Limit 

1955  OUTPUT  KBD; Cl  ear® ; 

1 SG0  !  - 

1965  1  Xfold  =  X  direction  e-folding  scale 

1970  PRINT  TABXY( Tx . Ty ), "Enter  the  X  direction  E-folding  scale" 

1975  DISP  "Xfold?"; 

1980  IF  Lp_f lg  THEN  OUTPUT  KBD; Xfold; 

1985  INPUT  "" .Xfold 

1990  OUTPUT  KBD; Clear®; 

1995  ! 

Z000  I  - 

Z005  1  Yfold  =  Y  direction  e-folding  scale 

Z010  PRINT  TABXY( Tx , Ty )," Enter  the  Y  direction  E-folding  scale.” 

Z015  DISP  "Yfold?"; 

Z0Z0  IF  Lp_f lg  THEN  OUTPUT  KBD; Yfold; 

Z0Z5  INPUT  "".Yfold 

Z030  OUTPUT  KBD; Clear®; 

Z035  ! 

Z040  I  - 

Z045  1  Tfold  =  Time  e-folding  scale 

Z050  PRINT  TABXY( Tx , Ty ), "Enter  the  Time  E-folding  scale  (seconds)." 

Z055  DISP  "Tfold?"; 

Z060  IF  Lp_f lg  THEN  OUTPUT  KBD; Tfold; 

Z065  INPUT  "".Tfold 

Z070  OUTPUT  KBD; Clear®; 

Z075  i  - 

Z080  !  Data_file$  =  Observed  data  file. 

Z085  PRINT  TABXY( Tx ,Ty) ."Enter  the  observed  data  file  including  file  specifier 
s . 

Z090  DISP  "Data_f ile$?" ; 

Z095  IF  Lp_f lg  THEN  OUTPUT  KBD; Dat a_f i 1 e$ ; 

Z100  LINPUT  "",Data_f ile® 

Z 105  Data_f ile$=UPC$( TRIM$( Data_f ile$) > 

Z 1 10  OUTPUT  KBD; Clear®; 

Z  1  1  5  !  - 

Z 1 Z0  !  Ip_file®  =  Interpolation  position  data  file. 

Z1Z5  PRINT  TABXY( Tx ,Ty) ."Enter  the  interpolation  position  file  including  file 
specif iers . " 

Z 1 30  DISP  "Ip_f ile$?"  ; 

Z 1 35  IF  Lp_f lg  THEN  OUTPUT  KBD; I p_f 1 1 e® ; 


2140  LINPUT  ,,M  ,  Ip_f  lie® 

2  145  Ip_f ile$=UPC$< TRI M$( Ip_f ile$>  > 

2150  OUTPUT  KBD; Cl ear$ ; 

2155  i  - 

2  160  1  Soa_fcst$  =  SOA  forecast  field  file. 

2165  PRINT  TABXY( Tx , Ty ) , " E nt er  the  file  specifier  for  OR  forecast  fields." 

2  170  DISP  " Soa_f  cst  $7 " ; 

2  175  IF  Lp_f lg  THEN  OUTPUT  KBD; Soa_f cst$ ; 

2180  LINPUT  "“,Soa_fcst$ 

2 1 85  Soa_f  cst$=UPC$( TRI M$(  Soa_f cst$> ) 

2190  OUTPUT  KBD; Clear? ; 

2  195  I  - 

2200  1  Soa_evar$  =  SOA  expected  error  file. 

2205  PRINT  TABXY( Tx , Ty ) , " E nt er  the  file  specifier  for  OR  expected  error  fields 


22  10 
2215 
2220 
2225 
2230 
2235 
2240 
2245 
2250 
2255 
2260 
2265 
»  #  " 

2270 

* 

2275 
2280 
2285 
2290 
2295 
2  300 
2  305 
2310 
2315 
2320 
2325 
2  330 
2335 
Z  340 
2  345 
2  350 
Z  355 
2  360 
Z  365 
2  370 


DISP  "Soa_evar$7" ; 

IF  Lp_f lg  THEN  OUTPUT  KBD; Soa_evar$; 

LINPUT  "",Soa_evar$ 

Soa_evar$=UPC$<  TRIM$(  Soa_evar$ ) ) 

OUTPUT  KBD; Cl ear$ ; 

I  - 

! 

1  -  Echo  check 

! 

T  ab=  1  0 

PRINT  TABXY( 10,2 ) 

PRINT  TAB( Tab) , "***********»  Echo  check  of  input  variables  #*»*♦***#* 
PRINT  TRB(  Tab ),  "Tine  of  the  analysis:  ";Date0$;"  ";Time0 


PRINT  TRB( Tab) , "Number  of  objective  analysis  to  make: 
PRINT  TAB< Tab ) , " T ime  interval: 

PRINT  TRB( Tab ) , "Max .  distance  radius: 

PRINT  TAB( Tab ) , "Max .  tine  radius: 

PRINT  TAB< Tab ), "Max inun  tine  lag: 

PRINT  TAB( Tab ), "Max inun  space  lag: 

PRINT  TAB( Tab ), “Max inun  number  of  influential  points: 
PRINT  TAB.(Tab),"X  direction  E-folding  scale 
PRINT  TAB(Tab),"Y  direction  E-folding  scale 
PRINT  TAB( Tab) ."Tine  E-folding  scale 
PRINT  TAB( Tab) ,"The  observed  data  file: 

PRINT  TAB( Tab ) , " The  interpolated  positions  file: 

PRINT  TRB( Tab ) , " The  SOR  forecast  output  file: 

PRINT  TRB( Tab ) , " The  SOR  expected  error  output  file: 
Lp_f lg=i 

Ansu$=" " 


" ;Nobj 
" ; Delta_t 
"  ;  XI  mit 
" ; T1 init 
" ; Max_t ime_lag 
" ; Max_space_l ag 
" ; L init 
" ; Xf old 
" ; Yf old 
" ; Tf  ol d 
" ; Data_f i 1 e$ 

" ; lp_f i 1 e$ 

" ; Soa_f cst$ 

" ; Soa_evar$ 


INPUT  "Do  you  want  to  change  an«  values  (Y/N)  <no>",Ansui$ 
EXIT  IF  UPC$<  RnsuSM  ;  1  ]  )<>"Y" 

END  LOOP 


I 
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2-575  OUTPUT  KBO;Clear$; 

2  380  i 

2385  1  read  in  the  observed  data 

2  390  i 

2  395  Get  _dat  a(  Dat a_f  i  I  e$  .  Xdat  a<  *  )  ,  Ydata<  *  )  .  T data!  *  )  ,  IJhidat  a(  *  )  ,  Ndata ) 

2400  IF  Ndata=0  THEN 

2405  PRINT  TA8( Tab )." Error  in  reading  the  observed  data;  Program  stopped" 

2410  STOP 
2415  END  IF 

2420  PRINT  TAB( Tab ), "Total  number  of  the  observed  data:  ".Ndata 
2425  ! 

2430  1  read  in  the  interpolated  positions 

2435  i 

2440  Get_ip( I p_f i le$ , Xip< * ) ,  Y  ip( * ) .Nip) 

2445  i 

2450  IF  Nip=0  THEN 

2455  PRINT  TAB( Tab) . "Error  in  reading  interpolated  positions;  Program  stopped" 
24G0  STOP 
2455  END  IF 

2470  PRINT  "Total  number  of  the  interpolated  pos it  ions : " , Nip 
2475  PRINT  "Open  " , Soa_f cst$ , “  for  forecast  output." 

2480  ASSIGN  @Fcst_out  TO  Soa_fcst$ 

2485  PRINT  "Open  " , Soa_evar$ . "  for  expected  error  output." 

2490  ASSIGN  @Evar_out  TO  Soa_evar$ 

2495  WAIT  1 
2500  i 

2505  T0=DATE(Date0$)+TIME(Time0S)  ‘  Convert  the  real  time  to  HP  time  format 
2510  ! 

2515  1  Loop  for  each  analysis 

2520  i 

25Z5  FOR  I obj  =  1  TO  Nobj 
2530  OUTPUT  KBD; Clear®; 

2535  Tf  =  T0+Del ta_t  *  I  obj 

2540  Tl=Tf-Tl imit 

2545  Tu=Tf  +  T1 imit 

2550  IF  Del ta_t=0  THEN 

2555  PRINT  TABXYC Tx , Ty ) , " Time  of  the  analysis  to  make:  " ,DATE$( Tf  > ,TIME$( Tf ) 
2560  ELSE 

Z5G5  PRINT  TABXY( Tx ,Ty) ."Forecast  time:  " ,DATE$( Tf  > ,TIME$<  Tf  > 

Z570  END  IF 

2575  ! 

Z580  !  get  the  observed  data  points  within  the  limited  range  and  at  the  proper 

time 

2585  ' 

2590  Get_rd( Xdata(  ») , Ydata( »  > ,Tdata(*),Phidata<*), Ndat a . T1 ,Tu,Xlimit,Xrd(*).Yr 
d( * ) . Trd( * ) .Phird( * ) ,N) 

2595  i 

2G00  Evar=0. 

2G05  IF  Evar> 1 .0  THEN 

2610  PRINT  TAB< Tab ) , "The  error  noise  level  100X  exceeded" 
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2615  STOP 

2620  END  IF 
Z625  ' 

2630  1  Do  the  analysis  for  each  point 

ZG35  ' 

2640  FOR  Ip-1  TO  Nip 

ZG45  COLL  Select! Limit, Xip( Ip), Yip( Ip) ,Tf,Xrd(*>,Yrd(*),Trd(*),Phird<*),Xopd( 

* ) , Yopd( *  )  ,Topd( * ) ,  Phiopd!  * ) , N , Nobs , Max_t ime_l ag , l“tax_space_lag ) 

Z650  IF  Nobs<>0  THEN  2670 

ZG55  Theta! Ip  >  =  99999 

ZG60  Es( Ip)  =39999 

ZGG5  GOTO  Next_ip 

2670  CALL.  Seal  ar_oa(  Xopd<  *  )  ,Yopd(  *  )  ,Topd(  *  >  ,Phiopd(  *  )  ,Nobs,X,Y,Tf  , Valp , U , I er  ) 

ZG75  ! 

2680  !  if  the  error  field  is  too  large,  increment  error  variance  and  re_do 

2685  i 

Z690  IF  I er=0  THEN  GOTO  2705 

Z695  E var=Evar+ . 0 1 

2700  GOTO  2605 

2705  Theta! Ip)=Valp 

Z710  Es(  Ip  >=SQR(  ABS!  Ul )  ) 

21 1 5  DISP  USING  "K , ZZZZ , K ,  MO. DDDE , K ,MD. DDDE "Point  no.  ".Ip,"  Theta=" , Theta 

(Ip),"  Expected  Errors-" , Es( I p ) 

2720  Next_ip : NEXT  Ip 
2725  i 

2730  Tab- 10 

2735  PRINT  TAB! Tab) 

2740  PRINT  TREK  Tab) , "The  diagnostics  of  the  observed  fields" 

2745  CALL  Diag( Phird< * ) , N ) 

2750  PRINT  TAB!  Tab ) , 

2755  PRINT  TAEK  Tab ) , "the  diagnostics  of  the  predicted  fields" 

2760  CALL  Diag<  Theta( * ) ,Nip> 

2765  PRINT  TAEK  Tab ) , 

2770  PRINT  TAEK  Tab ) , " the  diagnostics  of  the  error  fields" 

2775  CALL  Di ag( Es( * ) , Nip ) 

Z780  I 

2785  1  output  the  O.A  .field 

2790  i 
Z795  BEEP 
2  800  ! 

2805  1  -  output  the  0.  A.  forcast  field. 

2810  i 

2815  DISP  "Output  the  Forecast  field  to  ",Soa_fcst$ 

2820  OUTPUT  @Fcst_out ; "0.  A.  Forecast  Field" 

2825  OUTPUT  @Fcst_out;"Time:  "&DATE$<  Tf  )8>"  "&TIME$(  Tf  ) 

2830  FOR  I p= 1  TO  Nip 

Z835  OUTPUT  @Fcst_out ; Thetaf Ip) 

2840  NEXT  Ip 
2845  ' 

2850  1  -  output  the  Expect  Error  f leld. 
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Z855  i 

Z860  DISP  "Output  the  error  field  to  " ,Soa_evar$ 

Z8G5  OUTPUT  @Evar_out ; "Expected  Error  Field" 

Z 870  OUTPUT  @Evar_out; "Time:  "&DATES! Tf >&"  "&TIMES! Tf ) 

Z 875  FOR  I p= 1  TO  Nip 

Z880  OUTPUT  @Evar_out;Es! Ip) 

Z885  NEXT  Ip 
Z890  DISP 
Z895  ' 

Z900  NEXT  Iobj 
2905  I 

Z9t0  1  close  all  output  files. 

Z915  ' 

Z9Z0  ASSIGN  @Fcst_out  TO  * 

Z925  ASSIGN  @Evar_out  TO  * 

Z930  Finish:  1 
Z935  DISP  "Finished" 

Z940  END 
Z945  ! 

2950  ' 

Z955  i 

Z9G0  Get_rd: SUB  Get_rd<  Xdata<  * ) , Ydata( * ) ,Tdata<  *  > ,Phidata( *>,M,Tl,Tu,Xl imit ,Xrd( 
* ) , Yrd< *  )  , Trd( * ) ,Phird<  * ) , N ) 

Z9G5  I 

Z970  1  routine  gets  the  data  before  and  after  a  given  time  and  uithin 

Z975  !  a  given  spatial  radius  from  the  domain  reference  point. 

Z980  i 
Z985  RAD 
Z990  OPTION  BASE  1 

Z995  IF  Tl=Tu  THEN 

3000  PRINT  USING  "K  .  K  ,  K“ ! "Use  data  on  "  ,DATE$<  Tu  >  ,  "  ",TIME$!Tu> 

3005  ELSE 

3010  PRINT  USING  "K,K,K,K,K,K,K,K";"Use  data  after  " , DATE$( T1 ) , "  TIMES! T1 ) , 
"  and  before  ",  DATES! Tu ), “  ",TIME$!Tu) 

3015  END  IF 

30Z0  PRINT  USING  " K , MD. DDDE , K " ; "  and  within  the  range  of  ".Xlimit,"  from  the  i 
ef erence  point. " 

30Z5  N=0 

3030  FOR  1=1  TO  M 

3035  IF  Tdata! I ) >Tu  THEN  GOTO  Next_i 

3040  IF  Tdata!  I  XT1  THEN  GOTO  Next_i 

3045  IF  ABS! Xdata! I ) )>Xlimit  THEN  GOTO  Next_i 

3050  IF  ABS! Ydata! I ) )>Xlimit  THEN  GOTO  Next_i 

3055  N=N+ 1 

30G0  Xrd! N )=Xdata< I ) 

3065  Yrd! N )=Ydata<  I  ) 

3070  Trd! N)=Tdata(  I  ) 

3075  Phird! N ) =Phidata(  I  > 

3080  Next_i : NEXT  I 
3085  IF  N=0  THEN 


B-33 


3090  PRINT  “No  data  found;  Please  re-define  tine  or  distance  radius;  Program 
stopped . " 

3095  BEEP 

3100  STOP 

3105  ELSE 

3110  PRINT  "Number  of  data  uithin  this  range:", N 

3115  END  IF 

3120  SUBENO 
31 Z5  * 

3130  ' 

3135  ! 

3 1 40  Select: SUB  Select! Limit ,X,Y , Tf ,Xrd( * ) ,Yrd( *  > ,Trd( * ) ,Phird( * ) ,Xopd( *  )  ,  Yopd( * 
) ,Topd( *  )  ,Phiopd( * ) ,N, Nobs . Max_t ime_l ag , Max_space_l ag ) 

3145  * 

3150  1  -  routne  eliminates  the  distant  (in  space  and  time)  points  and  selects 

3155  !  the  most  "LIMIT"  near  points  to  an  interpolation  point  X.Y.Tf. 

3160  I 

31 B5  OPTION  BOSE  1 
3170  ROD 

3175  DIM  Index( 1089 ) ,Cor<  1089 ) 

3180  COM  /Fold/  Xfold.Yfold.Tfold 

3185  COM  /Efield/  Evar 

3190  COM  /Cblock/  C(40) 

3195  Cphase=0. 

3200  Nobs=0 

3205  FOR  1=1  TO  N 

3210  Delt=Tf -Trd< I ) 

3215  IF  ABS( Del t ) >Max_t ime__l ag  THEN  GOTO  Next_i 

3220  Delx=X-Xrd( I ) 

3225  Oel y=Y'Yrd( I ) 

32  30  R=SQR( ( Delx ~Cphase*Del t ) * ( Del x-Cphase*Del t )  +  Del y*Del y  > 

3235  IF  R>Max_space_lag  THEN  GOTO  Next_i 
3240  Nobs=Nobs+1 

3245  I ndex( Nobs )=I 

3250  Cor( Nobs  >=FNCov( Delx ,Dely,Delt , Xfold.Yfold.Tfold) 

3255  Next_i : NEXT  I 
3?S0  IF  Nnbs=0  THEN 
32B5  BEEP 

3270  OUTPUT  KBD; CHR$<  255  ) ; "K"  ; 

3275  PRINT  TABXY( 10, 10) ,"***  Warning  .^rningl  Warning!  ***" 

3280  PRINT  USING  "K,K,K";"No  data  were  selected  for  interpolated  position:  ", 

X ,  Y 

3285  PRINT  "The  present  radii  of  the  influence  domain  are:" 

3290  PRINT  "Maximum  space  lags:  " , Max_space_l ag 

3295  PRINT  "maximum  time  lags:  " , Max_t ime_l ag 

3300  PRIN1  'User  responses:" 

3305  PRINT  "  1)  Use  larger  maximum  time  or  space  lags,  then  start  over." 

3310  PRINT  "  2)  Assign  a  value  of  99999  to  the  estimated  value  at  this  posi 

tion  and  go  to  the  next  point." 

3315  INPUT  "Enter  your  option  (1  or  2)  <2>", Option 
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3320  OUTPUT  KBD; CHR$( 255 ) ; "K" ; 

3325  IF  Opt  l on= 1  THEN  STOP 

3330  SUBEXIT 

3335  END  IF 

3340  IF  Nobs>Limit  THEN 

3345  1  REDIM  Cor( Nobs  )  , I ndex ( Nobs  ) 

3350  '  MAT  SORT  Cor< * >  DES  TO  Index 
3355  Sort ( Cor (*), I ndex (*), Nobs ) 

3350  Nobs=Limit 

3355  END  IF 

3370  FOR  1=1  TO  Nobs 

3375  J=I ndex( I  ) 

3380  Xopd< I )=Xrd( J ) 

3385  Yopd( I  )=Yrd< J) 

3390  Topd< I )=Trd<  J ) 

3395  PhiopdT I )=Phird< J) 

3400  C( I >=Cor( I  ) 

3405  NEXT  I 
3410  SUBEND 
3415  1 

3420  ! 

3425  i 

3430  Scalar_oa: SUB  Scalar_oa<  Xopd<  * ) , Yopd( *  > ,Topd( * ) ,Phiopd( * ) ,N 
er  ) 

3435  i 

3440  !  -  The  scalar  space-time  objective  analysis  routine. 

3445  i 

3450  OPTION  BASE  1 

3455  RAD 

3450  DIM  A< 40,40) 

3465  COM  fC block/  C<40) 

3470  COM  /Efield/  Eva r 

3475  IF  N<=0  THEN  GOTO  3BZ5 

3480  CALL  Set_i nva( A( * ) , Xopd( * > , Yopd( * ) , Topd( * > , N , I er  ) 

3485  IF  I er >0  THEN  SUBEXIT 

3490  CALL  Est_mean(  A< * ) ,Phicpd< * ) ,N, Ave ,SumZ > 

3495  i 

3500  I  U  is  the  error  variance 
3505  ! 

3510  W=0. 

3515  WZ=0. 

35Z0  I 

35Z5  1  -  valp  is  the  inpterpol ated  data 

3530  ' 

3535  Valp=0. 

3540  FOR  1=1  TO  N 

3545  H=0. 

3550  Dumc=C( I ) 

3555  FOR  J=1  TO  N 

3560  P=Dumc*C( J)*A< I , J) 


X  ,  Y  ,  Tf  ,  Valp  .  LI ,  I 
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3565  W=W+P 

3570  PZ  =  A< I . J)»Phiopd( J  ) 

3575  H=H+PZ 

3580  P3=C( J  >  *A( J , I) 

3585  U  Z  =  U  Z  *■  P  3 

3590  NEXT  .J 

3595  Oumy=Oumc*H 

3600  Val p= Val p+Dumy 

3605  NEXT  I 
3510  Valp=Vaip+Ave 
3615  Um=( 1 . - W2 ) A Z/ Sum2 
36Z0  W=( 1 . -  W ) +  Um 

36Z5  SUBEND 
3630  i 
3635  i 
3640  ! 

3645  Set_iriva:SU8  Set_inva(A(*),Xopd(*),Yopd(*>,T  opd(  *  )  ,  Nobs  ,  I  er  > 

3650  i 

3655  I  -  routine  sets  up  the  correlation  function  for  the  observations 

3660  •  given  the  positions  Xopd,  Yopd  and  tines,  topd,  it  returns  the 

3665  I  inverted  correlation  function  matrix. 

3670  ! 

3675  OPTION  BASE  1 

3680  RAD 

3635  DIN  Ip( 1089 ) 

3690  COM  /Fold/  Xf old , Yf old , Tf old 

3695  COM  /Efield/  Ever 

3700  Guard= 1 . 0E-30 

3705  ler=0 

3710  FOR  1=1  TO  Nobs 

3715  FOR  J=I  TO  Nobs 

37Z0  Delt=Topd( I )-Topd< J) 

37Z  5  Oelx  =  Xopd( I )-Xopd( J) 

3730  Dely=Yopd( I )-Yopd( J > 

3735  A<  I  ,  J)=FNCov(  Del  x  ,  Del y  ,  Del  t  ,  Xf  ol d  ,  Yf  ol  d  ,  Tf  ol  d ) 

3740  A( J , I ) =A( I , J ) 

3745  NEXT  J 

3750  A< I , I >  =  A< I ,1 )+Evar 

3755  NEXT  I 

3760  I 

3765  1  invert  the  matrix  a 

3770  i 

3775  CALL  lnvmtx( A( * > ,Nobs , A( » ) ,Nobs ,Nohs ,D, Ip(  * ) , Ier ) 

3780  IF  0<Guard  THEN 

3785  PRINT  "***  WARNING  THE  DETERMINANT  IS  VERY  SMALL;  DET=" , D 

3790  PRINT  "Suggestion  to  user:  Use  smaller  max.  influential  points." 

3795  I  er  ■=  -  1 

3800  END  IF 

3805  SUBEND 

3810  i 
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3815 

i 

3820 

i 

38Z5 

Sort : SUB  Sort( Cor( *) , Index! *) ,  N  ) 

3830 

1 

3835 

i  -  routine  sorts  the  index  and  correlation  in 

descending  or 

3840 

1 

3845 

OPTION  BASE  1 

3850 

RAD 

3855 

Igap=N 

3850 

IF  Igap<= 1  THEN  SUBEXIT 

3855 

Igap=Igap/2 

3870 

I  Flax  =N- I  gap 

3875 

I  ex  =0 

3880 

FOR  1=1  TO  I Fiax 

3885 

Iplusg=I+Igap 

3890 

IF  Cor< I ) >=Cor( Iplusg )  THEN  GOTO  Next_i 

3895 

Save=Cor< I ) 

3900 

Cor! I )=Cor( Iplusg) 

3905 

Cor! Iplusg)=Save 

39  10 

I save= I ndex ( I ) 

3915 

I ndex ( I )  =  I ndex  < Iplusg) 

3920 

Index! Ipl usg )  =  1  save 

3925 

I  ex=  1 

39  30 

Next_i : NEXT  I 

39  35 

IF  I ex <>0  THEN  GOTO  3875 

3940 

GOTO  3860 

3945 

SUBEND 

3950 

3955 

3960 

3965 

! 

I 

i 

Est_Fiean:  SUB  Est_Fiean(  A(*),Psi(*),N,Ave,  SumZ  ) 

3970 

3975 

i 

! - routine  calculates  the  estimated  Fiean,  and 

then  removes 

3980 

!  the  estimated  mean  from  the  observations  array. 

3985 

i 

3990 

OPTION  BASE  1 

3995 

RAD 

4000 

DIM  C< 1 089 ) , D( 1089) 

4005 

FOR  1=1  TO  N 

4010 

C< I >=0. 

4015 

D(  I )=0. 

4020 

FOR  K= 1  TO  N 

4025 

C( I )=C< I )+A( I ,K )*Psi( K ) 

4030 

D( I )=D< I >+A! I ,K  > 

4035 

NEXT  K 

4040 

NEXT  I 

4045 

Sum  1 =0. 

4050 

Sum2=0. 

4055 

FOR  1=1  TO  N 

4060 

Sum1=Sum1 +C( I  ) 
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4065  SumZ  =  SumZ  iD(  1  ) 

4070  NEXT  I 
4075  i 

4080  1  -  calculate  the  estimated  mean 

4085  i 

4090  Ave=Sum 1 /SumZ 
4095  i 

4100  !  -  remove  the  estimated  mean  from  the  observations. 

4105  I 

4110  FOR  1=1  TO  N 
4115  Ps i ( I ) =Ps i ( I ) - Ave 
41Z0  NEXT  I 
41Z5  SUBEND 
4130  ' 

4135  l 
4140  ! 

4145  InvmtxiSUB  I nvmt x < A( * ) , Na , * > , Nv , N , D , I p( * ) , I er ) 

4150  ! 

4155  >  -  routne  inverts  the  Matrix  A. 

41G0  1 

4 1 G5  1  V  is  the  inverted  matrix  of  A. 

4170  !  D  is  the  determinent  of  Matrix  a. 

4175  i 

4180  OPTION  BASE  1 

4185  RAD 

4190  Iexmax=75 

4195  Ier  =  FNIerinv<  N.Na.Nv) 

4Z00  IF  Ier<>0  THEN  4555 

4Z05  FOR  J=1  TO  N 
4Z10  Ip(J)=0 

4Z1S  FOR  1=1  TO  N 

4ZZ0  V< I , J)=A( I , J) 

4ZZ5  NEXT  I 

4Z30  NEXT  J 

4Z35  D=  1  . 

4Z40  lex=0 

4Z45  FOR  M= 1  TO  N 
4Z50  Vmax=0. 

4Z55  FOR  J=1  TO  N 

4ZG0  IF  Ip(J)<>0  THEN  4305 

4ZG5  FOR  1=1  TO  N 

4Z70  IF  Ip<I><>0  THEN  4300 

4Z75  Vh=ABS<  V( I , J ) ) 

4Z80  IF  Vmax^Vh  THEN  4300 

4Z85  Vmax=Vh 

4Z90  K=I 

4Z95  L= J 

4300  NEXT  I 

4305  NEXT  J 

4310  I p<  L ) =K 
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4315  Npm=N+M 

43Z0  I p<  Npm ) =L 

43Z5  D=D* V( K . L  > 

4330  IF  ABS(D)  =1.0  THEN  4350 

4335  D=D* . 1 

4340  Iex=Iex+1 

4345  GOTO  4330 

4350  Pvt=V(K,L> 

4355  IF  M=  t  THEN  Pv tmx=ABS<  P vt ) 

4360  IF  ( ABS( Pvt/M)+Pvtmx )=Pvtmx  THEN  45Z0 

4365  V(  K  , L  )  =  1  . 

4370  FOR  J=1  TO  N 

4375  Hoi d=V( K  ,  J ) 

4380  V<K,J)=V<L.J> 

4385  V( L , J  )  =  Hol d/Pvt 

4390  NEXT  J 

4395  FOR  1=1  TO  N 

4400  IF  I =L  THEN  4430 

4405  Hold=V( I ,L  > 

4410  V( I , L ) =0. 

4415  FOR  J=1  TO  N 

44Z0  V  < I ,J)=V(IfJ)-V(L,J)*Hold 

44Z5  NEXT  J 

4430  NEXT  I 

4435  NEXT  M 

4440  M=N+N+ 1 

4445  FOR  J=1  TO  N 

4450  M=M- 1 

4455  L=Ip( M ) 

4460  K=Ip( L ) 

4465  IF  K=L  THEN  4500 

4470  D=-D 

4475  FOR  1=1  TO  N 

4480  Hoi d=V( I , L ) 

4485  V( I  ,L)=V( I  ,K) 

4490  V< I ,K  )=Hold 

4495  NEXT  I 

4500  NEXT  J 

4505  IF  I  ex >1  exmax  THEN  4540 
4510  D=D* 10. *  lex 
4515  GOTO  4555 
45Z0  Ier=33 
45Z5  BEEP 

4530  PRINT  TflB( 1 5  > , "UftRNING:  MATRIX  SINGULAR  IN  INVMTX" 

4535  GOTO  4555 

4540  Ier= 1 
4545  D=Iex 

4550  PRINT  TAB( 15) ."WARNING:  DETERMINANT  TOO  LARGE  IN  INVMTX" 
4555  SUBEND 
4560  l 


B-39 


4565  Fruerinv:DEF  FNI  er  mv(  N  ,  Na  ,  Nv  ) 

4570  OPTION  BASE  I 

4575  RAD 

4580  lerinv=0 

4585  IF  N>=  t  THEN  4605 

4590  Iennv  =  34 

4595  PRINT  "N  <  I  IN  INVMTX" 

4600  GOTO  4640 

4605  IF  Na>=N  THEN  4625 

4610  Iennv=35 

4615  PRINT  "NA<N  IN  INVMTX" 

4620  GOTO  4640 

4625  IF  Nv>=N  THEN  4640 

4630  Ierinv=36 

4635  PRINT  "NV  <  N  IN  INVMTX" 

4640  RETURN  Ierinv 

4645  FNEND 
4650  ' 

4655  ! 

4660  ! 

4665  DiagtSUB  Oiag( Dps i( * > , M ) 

4670  i 

4675  !  routine  calculates  the  statistical  parameters  of  input  array 

4680  !  such  as  the  mean,  variance,  root  mean  square,  mimimum  and 

4685  !  maximum. 

4690  ! 

4695  OPTION  BASE  1 

4700  RAD 

4705  ! 

4710  IF  M<=0  THEN  SUBEXIT 

4715  PRINT  "Number  of  points:  ",M 

4720  Ave=0. 

4725  Sdv=0. 

4730  Psimax  =  Dpsi(  1  ) 

4735  Psimin=Dpsi< 1 ) 

4740  Ave=Dps i  (  1  ) 

4745  IF  M=1  THEN  Cal_rms 

4750  FOR  1=2  TO  M 

4755  IF  Psimax<Dpsi(  I  )  THEN  Psimax-.'  i<  I  > 

4760  IF  Psimin>Dpsi( I )  THEN  Psimin=Dpsi< I > 

4765  Sdv=( < 1-2 )«Sdv+< 1-1 )*( Dpsi( I )-Aver ) "2/ I >/< I- 1 ) 

4770  Ave=( ( I -  1 ) * Ave+Dpsi( I ) ) / I 
4775  NEXT  I 

4780  Cal_rms : Rm5=SQR(  Sdv+Ave',2  > 

4785  PRINT  "Mean  and  Variance:  ",Ave,5dv 
4790  Sdv=SQR<  Sdv ) 

4795  PRINT  "Standard  Deviation:  ",Sdv 

4800  PRINT  "RMS  of  field:  ",Rms 

4805  PRINT  "Minimum  of  field:  ".Psimin 

4810  PRINT  "Maximum  of  field:  ", Psimax 
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43  IS  SUBEND 
4820  i 
4825  ' 

4830  I 

4835  be  t  _  dat a:  SUB  Get  dat a( F i 1 enameS .  Xdata(  *  )  .Ydata(  *  >  ,Tdata(  *  )  ,Phi dat a(  »  )  , Ndat  a 

) 

4840  I 

4845  1  -  routine  reads  in  the  observed  data. 

4850  I 

4855  ASSIGN  @Path_in  TO  Filename® 

4860  Ndat  a=0 

4865  PRINT  "  »*»**«•*»**  Echo  check  of  the  first  ten  records  »♦»*»»**«» 

*  " 

4870  PRINT  "  No.  X_P0S  Y_POS  TIME  DATA" 

4875  ON  END  @Path_in  GOTO  Close_file 

4880  DISP  USING  "K ,K ,K" ; "Reading  the  observed  data  from  ".Filenames,",  please 
wait.  " 

4885  LOOP 

4890  ENTER  @Path_in; X , Y , T ,Phi , Dummy 

4895  IF  Ndata<,=  1 0  THEN  PRINT  USING  "DDDO ,  XX ,  MO.  DODDDDE  ,  XX ,  MD.  DDDDDDE  ,  XX  ,  MD.  DO 

DDDDDDDE , XX , MD. DDDDDDE" ; Ndat a , X , Y , T , Phi 

4900  Ndata=Ndata+ I 

4905  Xdata( Ndata)=X 

4910  Ydata<  Ndata)=Y 

4915  Tdatal Ndata )=T 

49Z0  Phidata( Ndata)=Phi 

49Z5  END  LOOP 

4930  Close_f ile: ASSIGN  @Path  in  TO  *  !  Closing  input_file. 

4935  DISP 
4940  SUBEND 
4945  ! 

4950  I 
4955  i 

4960  Get_ip:SUB  Get_ip< F i 1 enameS , Xdat a< * > , Ydat a( * ) , Nip ) 

4965  l 

4970  1  -  routine  reads  in  the  computating  points. 

4975  ' 

4980  ASSIGN  @P<jth_in  TO  Filename® 

4985  PRINT  "*  Echo  check  of  the  first  ten  records  *" 

4990  PRINT  "  No.  X_IP_P0S.  Y_IP_P0S." 

4995  Nip=0 

5000  ON  END  @Path_in  GOTO  Close_file 

5005  DISP  USING  "K , K , K" ; "Reading  the  interpolated  positions  from  ".Filenames," 
,  please  wait." 

5010  LOOP 

5015  ENTER  @Path_in;X,Y 

50Z0  IF  N i p <  = 1 0  THEN  PRINT  USING  "DDDD , XX , MD . DDDDDDE , XX , MD. DDDDDDE” ; Nip , X , Y 
50Z5  Nip=Nip+1 

5030  Xdata( Nip  >  =  X 

5035  Ydatal Nip  >  =Y 


B-41 


5040  END  LOOP 

5045  C 1 ose_f l 1 e : ASSIGN  @Path_in  TO  *  1  Closing  input_file. 

5050  DI SP 
5055  SUBEND 

5060  OEF  FNCov<  Deix ,  Del  y ,  Del  t ,Xfold,Yfold,Tfold> 

50G5  i 

5070  1  routine  gets  the  correlation  function  from  a  fitted  formula. 

5075  i 

5080  OPTION  BOSE  I 
5085  ROD 

5090  COM  /Efield/  Evar 

5095  RETURN  EXP<  -  SQR<  (  Del x/Xf old > *( Oelx/Xf ol d >+( Del y/ Yf ol d ) * < Del y/ Yf ol d  )  +( Del t 
/Tfold)*( Delt/Tfold) ) ) 

5100  i  R2  =  SQR(  <Delx-'Z+0ely''2  )/5.0Et4AZ) 

5105  I  Cor=EXP( -RZ ) 

5110  I  RETURN  Cor 
5115  FNEND 
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1000 
1001 
1005 
1010 
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10Z0 
I  0Z5 
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1200 
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PROGRAM  VOA 

CHANGES  TO  PROGRAM  BY  A. ALLEN  AND  M. D. COUTURIER  16  DEC  1987 


*  » 

»  VECTOR  SPACE-TIME  OBJECTIVE  ANALYSIS  » 

*  * 

*  Language:  BASIC  3.0  * 

*  System:  Hewlett  Packard  9816  * 

*  Version:  1 . @0  * 

*  Date:  July  1986  * 

*  » 

*  Revisied  by  Dr.  L.  Charles  Sun  * 

*  Hawaii  Institute  of  Geophysics  * 

*  University  of  Hauaii  * 

*  Honolulu,  Hauaii  96822  * 

*  * 

*  Developed  at  U.  S.  Coast  Guard  Research  And  Development  Center  * 
»  Grot  in,  Connecticut  06340  * 

»  * 


i  ##»*#»*»*#***##»****##*****###**#*******#***#*#**»»**#»*»****«*»#*** 
i 

OPTION  BASE  1 
RAD 

DIM  Message®! 2561 

DIM  Xdata(  1089 ) , Ydata(  1089 ) ,Tdata<  1089  > ,Udata(  1089 ) ,Vdata( 1089) 

DIM  Xrd< 1089 ) ,Yrd< 1089) ,Trd< 1089 ) ,Urd< 1089) ,Vrd( 1089) 

□  I M  XopdC 1089), Yopd(  1089), Topd< 1089), Uopd< 1089), Vopd( 1 089 ) 

DIM  Xip< 1089) ,Yip< 1089) ,Tip( 1089) 

DIM  Uest< 1089 ) ,Vest< 1089) ,Uerror( 1089) ,Verror( 1089) ,Unrms< 1089) ,Vnrms( 1089 

COM  /Carr/  Cuu( 33 , 33 ) , Cuv< 33 , 33 ) , Cvv( 33 , 33 ) , Xbin , Ybin , Mbin , Nbin 
COM  /Fold/  Xfold, Yfold.Tfold 
COM  /Ef ield/  El ,EZ 
C 1 ear$=CHR$< 255)&"K" 

Sdate$=DATE$( TIMEDATE ) 

St ime$=TI ME$( TIMEDATE ) 

I 

OUTPUT  KBD?Clear$i  1  clear  the  screen. 

I 

PRINT  TABXY( 1 5 , 10) , "Do  you  need  documentations?" 

INPUT  "Enter  Y/N  for  yes/no  <no>",Ans$ 

OUTPUT  KBD;Clear$;  !  clear  the  screen. 

IF  UPC$< Ans$! 1 ; 1 1 >  <>" Y”  THEN  GOTO  14Z5 
Ty=  1  0 

PRINT  TABXY( 1 0 .7 ) , "Which  section  do  you  want  to  look  at  ?" 

PRINT  TAB( Ty) ."Section  1:  INTRODUCTION" 

PRINT  TAB(Ty),"  2:  BASIC  THEORY" 

PRINT  TAB( Ty ) , “  3:  DETERMINATION  OF  CORRELATON  FUNCTION" 

PRINT  TAB(Ty),"  4:  ELEMINATION  OF  DISTANT  DATA" 
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1 2  35  PRINT  TAB(Iy),"  5:  PROGRAMS  DESCRIPTIONS" 

1240  PRINT  TAEn  Ty  > , "  G:  USER  INSTRUCTIONS  AND  EXAMPLES" 

1245  PRINT  TAB(Ty),"  7:  REFERENCES" 

1250  PRINT  TAB(Ty),"  8:  DEFINITIONS  OF  VARIALBES" 

1255  PRINT 

1260  PRINT  TAB< Ty >," Enter  E  to  exit  help  nanus,  and  start  computation." 

1 2G5  INPUT  "Selection?" .Code® 

1Z70  OUTPUT  KBD; Clear®; 

1275  IF  Code$="F"  OR  Code$="e"  THEN  GOTO  1425 
1280  F i 1 e$="Doc_"&Code$ 

1285  ASSIGN  @Read_message  TO  File® 

1290  Read_message :  1 

1295  Count=0 

1 300  ON  END  @Read_message  GOTO  End_mess 
1305  ENTER  @Read_message; Message® 

1310  PRINT  TAB( 7 ) , Message®! 8 ; 801 
1315  Count=Count+ 1 
1320  IF  Count= 1 8  THEN 
t3Z5  INPUT  "More  message  ?  “,Ans® 

1330  IF  Ans$=" y"  OR  Ans$="Y"  THEN 
1335  OUTPUT  KBD; Cl  ear®; 

1 340  GOTO  Read_message 

1345  ELSE 

1350  OUTPUT  KBD; Clear®; 

1355  GOTO  1210 

1360  END  IF 
1365  END  IF 
1 370  GOTO  1 300 
1375  End_mess:  1 

1380  INPUT  "End  of  this  section,  Need  another  section?  ",Ans$ 

1385  IF  Ans$="y"  OR  Ans$="Y"  THEN 
1390  OUTPUT  KBD; Cl  ear®; 

1395  GOTO  1210 
1400  ELSE 

1405  OUTPUT  KBD; Clear®; 

1410  GOTO  1425 
1415  END  IF 
1420  i 

14Z5  I nput_parameter  :  ! 

1430  Tab=7 

1435  PRINT  TAB(Tab) 

1440  PRINT  TAB(Tab) 

1445  PRINT  TAB(Tab) 

1450  PRINT  TAB< Tab  I,"********************************************************** 
#**###*#*»*»” 

1455  PRINT  TAB< Tab) ,"* 

*  M 

1460  PRINT  TAB(Tab),"*  Vector  Space-Tine  Objective  Analysis  Package 

#  " 

1465  PRINT  TAB( Tab ) , " » 
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Language:  BASIC  3.01 


1470  PRINT  T AB(  Tab) , " » 

*  “ 

1475  PRINT  TAB(Tab),"*  System:  Hewlett  Packard  9000-2000 

#  " 

i480  PRINT  TAB(Tab),"*  Version  1.00 

*  " 

1485  PRINT  T  AB( Tab ) . " *  July  19  86 

#  " 

1490  PRINT  TAB( Tab) . “ * 

*  " 

1495  PRINT  TAB(Tab),"*  Coast  Guard  Research  and  Development  Center 

»  " 

1500  PRINT  TAB(Tab),"*  Groton,  Connecticut  06340 

#  " 

1505  PRINT  TAB< Tab) 

*  " 

1510  PRINT  TAB( Tab ),"***»*»***»*»****»******»*****************»****»**********» 
»****#*#***»" 

1515  UAIT  1 

1520  OUTPUT  KBD;Clear$; 

1525  PRINT 
1530  PRINT 
1535  Tab= 1 0 

1540  PRINT  TAB(  Tab ), "Before  you  go  on,  please  check  the  followings:" 

1545  PRINT  TAB(Tab),"  1)  Is  the  disc  containing  the  observed  data  and" 

1550  PRINT  TAB(Tab),"  the  inter/extraplotat ion  positions  in  DRIVE#  1 

1555  PRINT  TAB<Tab),“  2)  are  the  output  files  existed  7  ,  if  not,  create  them  b 

ef ore" 

1560  PRINT  TAB( Tab ) , "  you  do  the  analysis.” 

1565  PRINT  TAB(Tab),"  3)  The  file  length  for  each  output  should  be  2  records  1 
onger  than" 

1570  PRINT  TAB(Tab),"  the  total  inter/extrapolation  position  points." 

1575  PRINT  TAB(Tab),"  4)  The  current  version  allows  2000  input  data  and" 

1580  PRINT  TAB(Tab),"  1089  inter/extrapolat ion  points." 

1585  PRINT  TAB(Tab),"  5)  This  version  gets  the  correlation  function  from  " 

1590  PRINT  TAB( Tab ) , "  a  fitted  formuls." 

1595  PRINT  TAB( Tab) 

1600  PRINT  TAB( Tab )," Send  your  inquires  to" 

1605  PRINT  TAB<  Tab  > , "  Dr.  L.  Charles  Sun" 

1610  PRINT  TAB(Tab),"  Deptartment  of  Oceanography" 

1615  PRINT  TAB(Tab),"  University  of  Hawaii- 

1620  PRINT  TAB( Tab  > , "  Honolulu,  HI  96822" 

1625  PRINT  TAB< Tab) ,"  (808)348-7633" 

1630  INPUT  "Hit  (RETURN!  or  [ENTER]  to  continue" ,Answ$ 

1645  Tx= 1 0 

1650  Ty=  1  0 

1651  ! 

165Z  Lp_flg=0 
1653  LOOP 
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1654  OUTPUT  K BO; Clear®; 

!  656  '  - 

1660  1  Dat e0$ , T 1 me0$  =  time  of  the  analysis  to  make. 

IGG5  PRINT  TABXYI T* ,Ty) , "Enter  tine  of  the  analysis  to  make.IDD  MMM  Y Y , HH : MM : SS 
)" 

1670  PRINT  TflB(  Ty)  ."(enclosed  uith  quotations.)1' 

1671  DISP  " Dat e0$? , T l me0$? " ; 

I67Z  IF  Lp_flg  THEN  OUTPUT  KBD; . ;Date0$; . , . ;Tine0$; . ; 

1675  INPUT  *••'  .Date0$,Time0$ 

1680  Date0$=UPC$< Dats0$> 

1685  OUTPUT  KBD; Clear®, ■ 

1690  *  - 

1695  !  Nobj  -  number  of  objective  analysis  to  nake. 

1700  PRINT  TABXYt Tx , Ty ), "Enter  number  of  objective  analysis  to  nake." 

1701  DISP  "Nobj  ?” ; 

1702  IF  Lp_f lg  THEN  OUTPUT  KBD; Nobj ; 

1705  INPUT  "".Nobj 

1710  OUTPUT  KBD; Clear®; 

1715  i  - 

1720  1  Delta_t  =  tine  interval  for  each  extrapolation  in  tine. 

1725  PRINT  TABXY( Tx . Ty ), "Enter  tine  interval  for  extraploateion  in  tine" 

1730  PRINT  TAB(Ty),"  or  zero  for  instantaneoue  computation." 

1731  DISP  "Delta_t?" ; 

1732  IF  Lp_f lg  THEN  OUTPUT  KBD;Delta_t; 

1735  INPUT  "" ,Delta_t 

1740  OUTPUT  KBD; Cl  ear®; 

1745  i  - - - 

1750  1  XI imit  =  Max.  distance  radius  from  the  reference  point  of  the  domain. 

1755  PRINT  TABXY( Tx , Ty ), "Enter  Max.  distance  radius  from  the  reference  point  of 
the  domain." 

1760  PRINT  TAB(Ty),“All  data  uiithin  this  range  are  retained." 

1761  DISP  "Xlimit?"; 

1762  IF  Lp_f lg  THEN  OUTPUT  KBD;Xlimit; 

1765  INPUT  "" ,X1 imit 

1770  OUTPUT  KBD; Clear®; 

1775  i  - 

1780  !  Tlimit  =  Max.  time  radius  befor  and  after  the  time  of  the  analysis. 

1785  PRINT  7ABXY( Tx .Ty) ."Enter  the  max.  time  radius  befor  and  after  the  time  of 
the  analysis." 

1790  PRINT  TAB( Ty  > . "A1 1  data  within  this  range  are  retained." 

1791  DISP  "Tlimit?"; 

1792  IF  Lp_f lg  THEN  OUTPUT  KBD; T1 imit; 

1795  INPUT  "" ,T1 imit 

1800  OUTPUT  KBO; Clear®; 

1805  i  - 

1810  !  Max_space_l ag  =  maximum  space  lag. 

1815  1  Max_time_lag  =  maximum  time  lag. 

1820  PRINT  TABXY< Tx , Ty ) , "Set  the  influence  domain  for  computation  points," 

1825  PRINT  TAB(Ty),"the  maximum  time  lags  >=  the  multiplication  of  time  interva 
1" 
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18.30  T  T-TBi  T,  ■  .  "  arid  the  number  of  the  analysis  to  mak  e .  " 

1835  PRINT  iriB(  rv  i  . "  the  maximum  spae  lags  -=  the  sizes  of  the  computation  domm 

r, .  " 

1840  PRINT 

1345  PRINT  IrTBi  f  v  )  .  "Enter  the  max  imum  space  and  time  lags." 

184G  DI5P  " Mas _apace_ 1 ag? ,  Max_t ime  lag7" ; 

1847  IF  Lp  _f 1 g  THFN  OUTPUT  KBO; Max_space_ l ag ,Mag_t ime_lag; 

1350  INPU1  .  Ma-'.  _  space_  1  ag  ,  Max  t  ime_  l  ag 
1855  OUTPUT  KBD;Clear$; 

1  860  •  - - - * - -• - _  

1 8G5  i  Limit  -  Max.  number  of  influential  points. 

1870  PRINT  rOBXY(  I* ,  Ty  )  ,  "Enter  the  maximum  number  of  influential  points  " 

1871  DISP  "Limit?"; 

1872  IF  Lp  _f 1 g  THEN  OUTPUT  KBDiLimit; 

1875  INPUT  “".Limit 

1880  OUTPUT  KBD; Clear®; 

1881  i  - - 

1882  i  Xfold  =  X  direction  e-folding  (0.36)  scale  for  the  correlation  function 

1883  PRINT  FABXY( Tx , fy ), "Enter  the  X  direction  E-folding  scale" 

1884  DISP  "Xfold?"; 

1885  IF  Lp_f Ig  THEN  OUTPUT  K  80 ; X  f  o 1 d ; 

1 88G  INPUT  ""  .Xfold 

1887  OUTPUT  KBD; Clear®; 

,888  i  - - - 

1890  i  Yf old  =  Y  direction  e-folding  (0.3G)  scale  for  the  correlation  function 

1891  PRINT  TABXY( Tx , Ty ), "Enter  the  Y  direction  E-folding  scale" 

1892  DISP  " Yf old?" ; 

1893  IF  Lp_f 1 g  THEN  OUTPUT  KBO; Yf  old; 

1894  INPUT  ""  ,  Yfold 

1895  OUTPUT  KBD; Cl ear$; 

legs  l  - 

1897  I  Tfold  =  Time  (5EC0NDS)  e-folding  (0.3G)  scale  for  the  correlation  functi 
on 

1898  PRINT  TftBXYt Tx.Ty) ."Enter  the  TIME  e-folding  scale  (SECONDS)" 

1899  DISP  "Tfold?"; 

1900  IF  Lp_f lg  THEN  OUTPUT  KBD; Tfold; 

1901  INPUT  "” .Tfold 

1902  OUTPUT  KBD; Clear®; 

1 904  I  - 

1905  i  Data_file$,  Nskip  =  Observed  data  and  number  of  data  to  be  skipped. 

1906  PRINT  TABXY( Tx , Ty ). "Enter  the  observed  data  file  including  file  specifiers 


1907  DISP  " Dat a_f l 1 e$?" ; 

1908  IF  L p_ f 1 g  THEN  OUTPUT  KBD; Dat a_f 1 1 e$ ; 

1910  INPUT  ,Data_f ile$ 

1911  Dat a_f i 1 e$=UPC$<  Oat a_f i 1 e$ ) 

1912  OUTPUT  KBD; ClearS; 

1915  i  - 

1920  i  Ip_file$  =  Interpolation  position  data. 

1925  PRINT  TABXY( Tx , Ty ), "Enter  the  interpolation  position  file  including  file  s 


pec  if lers. " 

I  9  Z  6  DI5P  ■  Ip_  f  lie*’’"  ; 

1327  [F  LpJ'lq  THEN  OUTPUT  KBD ;  I  p_f  1 1  e$ ; 

19  30  INPUT . .  Ip_File® 

19  35  Tp_f iLe$-UPC$< Ip_f ile$> 

1940  OUTPUT  KBD; Clear®; 

1345  i  - 

1950  1  Voa_fcst$  =  VOA  forecast  field  file. 

1955  PRINT  TABXYI Tx , Ty ), "Enter  the  file  specifier  for  Oft  forecast  fields.” 

1956  DISP  " Voa_f cst$?" ; 

1957  IF  Lp_f Ig  THEN  OUTPUT  KBD; Voa. _f cst® ; 

I960  INPUT  "" , Voa_f cst® 

1 9G5  Voa_f cst $=UPC$(  Voa_f cst® ) 

1970  OUTPUT  KBD; Clear®; 

1975  '  - 

1980  1  Voa_evar$  =  VOA  expected  error  file. 

1985  PRINT  TABXY( Tx , Ty ), "Enter  the  file  specifier  for  Oft  expected  error  fields. 

1986  DISP  " Voa_evar$T " ; 

1987  IF  Lp_f 1 g  THEN  OUTPUT  KBD; Voa_evar$; 

1990  INPUT  "" ,Voa_evar$ 

1995  Voa_evar$=UPC$< Voa_evar $ ) 

Z000  OUTPUT  KBD; Clear®; 

Z005  i  - 

2010  i 

2015  1  -  Echo  check 

2020  l 

2025  T  ab= 1 0 

2030  PRINT  TABXY( 10,2 ) 

2035  PRINT  TAB< Tab ),"**#*#***#* **  Echo  check  of  input  variables  »***»***»»# 

*  " 

2040  PRINT  TAB< Tab )," Tine  of  the  analysis:  ",Date0$,"  ",Time0$ 

2045  PRINT  TfiB( Tab) , “Number  of  objective  analysis  to  make:  ",Nobj 

2050  PRINT  TAB( Tab) , " Tine  interval:  " ,Delta_t 

2055  PRINT  TAB( Tab ) , "Max .  distance  radius:  ",Xlimit 

20G0  PRINT  TAB( Tab ) , "Max .  tine  radius:  ".Tlinit 

2065  PRINT  TAB( Tab ) , "Max inun  tine  lag:  " , Max_t ine_l ag 

2070  PRINT  TAB< Tab ) , "Max inum  space  lag:  " , Max_space_l ag 

2075  PRINT  TftB( Tab ), "Max inum  number  of  influential  points:  ", Limit 

2080  PRINT  TflB<  Tab  )  ,  “The  observed  data  file:  "  ,Data_f  ile® 

Z085  PRINT  TAB( Tab ) , "The  interpolated  positions  file:  ",Ip_file$ 

2090  PRINT  TflB( Tab) ."The  VOA  forecast  output  file:  ",Voa_fcst$ 

2095  PRINT  TAB( Tab) , "The  VOA  expected  error  output  file:  ",Voa_evar$ 

2100  Lp_f lg= 1 

2101  flnsu$="" 

2102  INPUT  "Do  you  wnat  to  change  any  values  (Y/N)  <no>",Ansu® 

2103  EXIT  IF  UPC$< Ansu*[ 1 ; 1 1 ><>"Y" 

2104  END  LOOP 

21 15  OUTPUT  KBD; Clear®; 

2120  I 
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Z1ZS  1  input  the  observed  data 
2  t  30  ' 

Z  I  35'  Get  dat  a(  Oat  a  __f  i 1 e$ , Xdat a<  ♦  )  .  Ydat  a<  *  )  ,  Tdat  a(  *  )  ,  Udat  a  (  *  )  .  Vdat  a(  *  )  ,  Ndata  > 
Z140  IF  Ndata -0  THEN 

2  145  PRINT  "Error  in  reading  the  observed  data'1 
2150  STOP 
ZIS5  END  IF 

2  1  G0  PRINT  "Total  number  of  the  observed  data:  ".Ndata 
2  165  ' 

2170  i  read  the  interpolated  positions 
Z  175  I 

Z  t  80  OUTPUT  KBD; Clear®; 

Z 1 85  Get_ i p( Ip_f ile$,Xip< * ) , Yip<  * ) ,Nip> 

Z 190  i 

Z 1 95  IF  Nip=0  THEN 

ZZ00  PRINT  "Error  in  reading  interpolated  positon.  Program  stopped" 

ZZ05  STOP 
2Z10  END  IF 

ZZ15  PRINT  "Total  number  of  the  interpolated  pos it  ions : " , Nip 
2Z20  ! 

ZZZ5  T0=OATE( Date0$)+TIME< Time0$ )  1  Convert  the  real  time  to  HP  time  format 
ZZ30  ! 

2Z35  1  Loop  for  each  analysis 

Z  240  i 

Z 245  FOR  Iobj=1  TO  Nob] 

Z250  Tf =T0+Delta_t *Iobj 

ZZ55  T1 =Tf -T1 imit 

ZZ60  Tu=Tf  +  T1 imit 

ZZG5  PRINT  "Forecast  time:  ", DATE$( Tf ), TIMES! Tf  ) 

ZZ70  E I =0. 

ZZ75  EZ=0. 

ZZ80  ! 

ZZ85  1  get  the  observed  data  points  within  the  limited  range  and  at  the  proper 

time 

Z290  i 

ZZ95  CALL  Get_rd( Xdata( » ) , Ydata( * ) , Tdat a(  * ) , Udataf  * ) , Vdata(  *),Ndata,Tl,Tu,Xlim 
it ,Xrd( * ) ,  Yrd( * ) ,Trd( * ) ,Urd( * ) ,Vrd<  * ) ,N> 

Z300  PRINT  "Number  of  data  to  use  in  OA:  ",N 
Z  305  ! 

2310  CALL  Di ag<  Urd( * ) , N ) 

Z 3 1 5  CALL  Diag( Vrd( * ) ,N) 

17.17)  I 

Z3Z5  IF  El >1 .0  THEN 

Z330  °PINT  "The  error  noise  level  100%  exceeded;  ERROR  var iance=” ,E 1 
Z  335  STOP 
Z  340  .  ’  D  I F 

Z  345  • 

Z350  !  n  :he  analysis  for  each  point 

Z  35"  ' 

Z 360  FOR  I p= 1  TO  Nip 
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2365  CALL  Seleet(LimitlXip(Ip)!Yip(Ip>,Tf,Xrd(»),Yrd<»>,Trd<*),Urd(*),Vrd<*> 
Xopd<  *),Yopd(*).Topd(*).  Uopd(  *  )  ,  Vopd(  »  >  ,  N  ,  Nobs  ,  Max_space_l  ag  ,  Max_t  i  me_l  ag  ) 

2  370  CALL  Vect or_oa( Xopd( * ) , Y opd<  * ) , Topd( * ) , Uopd( * ) ,  Vopd( *  )  , Nobs ,Xip( Ip) ,Yip 
Ip) , Tf ,  Ues . Ves .Uerr.Verr.Ier) 

2  375  ' 

2380  IF  I er=0  THEN  2400 

2385  El =E I + . 01 

2390  E  2  =  E  2  +  .  0 1 

2395  GOTO  2325 

2400  Uest(Ip)=Ues 

2405  Vest ( I p ) =Ves 

2410  Uerror( I p ) =Uerr 

2415  Verrort Ip >=Verr 

2420  OISP  USING  " K , GOOD , K , MD. OOOE , K , MD. DDDE'V'Doi ng  V.O.A.  at  point  no.  ".Ip 

1  Est.  U:  "  ,Uest( Ip) Expected  Error : " , Uerrorf I p > 

2425  NEXT  Ip 
2430  i 

2435  !  compute  the  mean  and  variance  of  the  fields 

2440  ! 

2445  PRINT 

2450  PRINT  "the  diagnostics  of  the  forecast  fields" 

2455  CALL  Diag( Uest < * ) , Nip > 

2460  CALL  Diag< Vest< * ) , Nip ) 

2465  PRINT 

2470  PRINT  "the  diagnostics  of  the  error  fields" 

2475  CALL  Diag< Uerrorl * > , Nip) 

2480  CALL  Oiag< Verror< * ) , Nip ) 

2485  ! 

2490  !  output  the  V.O.A.  forcast  field. 

2495  i 

2500  OUTPUT  @Fcst_out ; "V.O. A.  forecast  field” 

2505  OUTPUT  @Fcst_out; "Time:  "&DATES< Tf >&"  "&TIME$<  Tf  ) 

2510  DISP  "Output  the  VOA  forecast  filed  to  ",Voa_fcst$ 

2515  FOR  I p= 1  TO  Nip 

2520  OUTPUT  @Fcst_outiUest< Ip) ,Vest< Ip > 

2525  NEXT  Ip 
2530  1 

2535  !  output  the  V.O.A.  error  variance  field. 

2540  i 

2545  OUTPUT  @Evar_out ; "Expected  Error  field" 

2550  OUTPUT  @Evar_out; "Time:  " &DATE$< Tf ) &"  "&TIME$( T f ) 

2555  DISP  "Output  the  VOA  error  variance  filed  to  ",Voa_evar$ 

2560  FOR  Ip=1  TO  Nip 

2565  OUTPUT  @Evar_out ; Uerror< I p ) . Ver ror< Ip ) 

2570  NEXT  Ip 
2575  i 

2580  DISP 
2585  NEXT  Iobj 
2590  ! 

2595  1  close  all  files 


B-50 


2600  ' 

2605  ASSIGN  @Fcst_out  TO  * 

Z6I0  ASSIGN  @Evar_out  TO  * 

2616  Cpu_t i me= TIME DATE- DATE ( Sdate$ ) - T I  ME ( Stime$> 

2620  PRINT  "CUP  time  for  this  run  =  " , Cpu_t ime 
2626  END 
2630  i 

2635  i - 

Z640  i 

2645  SUB  Get_rd(Xdata<  *),Ydata(*),T data! * ) , Udat a( * ) , Vdat a( *  ) ,M,T1 ,Tu,Xlimit ,Xrd 
( * ) , Yrd( *  > , Trd( * ) ,Urd( » ) ,Vrd< *  > ,N) 


2650 

I 

2655 

!  -  routine  gets  the  data  before  and  after  a  given  time  and  within 

2660 

!  a  given  spatial  radius  from  the  domain  reference  point. 

2665 

! 

2670 

RAD 

2675 

OPTION  BASE  1 

2680 

IF  Tl=Tu  THEN 

2685 

PRINT  USING  ”K  ,K  ,K” ; "Use  data  on  " , DATE$( Tu ) , " 

"  , TI ME$(  Tu) 

2690 

ELSE 

2695 

PRINT  USING  "K ,K ,K ,K  ,K  ,  K  ,K ,K" ;"Use  data  before 

" ,DATE$( Tu) " , TIME$( Tu ) 

and  after  " , DATE$( T1 ) , "  ",TIME$(T1> 

2700 

END  IF 

2705 

PRINT  USING  "K , MD. DDDE , K"  ;  "  and  within  thin  the 

range  of  ".Xlimit,"  from 

the  reference  point." 

2710 

N=0 

2715 

FOR  1=1  TO  M 

Z720 

IF  Tdatal I ) >Tu  THEN  GOTO  Next_i 

2725 

IF  TdatadXTl  THEN  GOTO  Next_i 

2730 

IF  ABS( Xdata( I ) )>Xlimit  THEN  Next_i 

2735 

IF  ABS(  Ydata<  I  )  )>Xlimit  THEN  Ne>.t_i 

2740 

N=N+ 1 

2745 

Xrd(N)  =  Xdata<  I  ) 

2750 

Yrd<  N)=Ydata< I ) 

2755 

Trd( N)=Tdata(  I ) 

2760 

Urd(N)=Udata(  I  ) 

2765 

Vrd( N)=Vdata<  I ) 

2770 

Next_i :  NEXT.  I 

2775 

IF  N=0  THEN 

Z780 

PRINT  "No  data  were  found  in  this  range,  please 

re-define  time  or  space 

ranges . " 

Z785 

PRINT  "Program  execution  halted." 

2790 

BEEP 

2795 

STOP 

2  800 

ELSE 

2  805 

PRINT  "No.  of  data  within  this  range:  ",N 

2810 

END  IF 

2815 

SUBEND 

2820  I 
2825  i 


B-51 


Z830  i 

Z835  SUB  Oe__mean(  Ph ird(  *  )  ,  M  ,  Mean  ,  Sdv  ) 

Z 840  OPTION  BASE  t 

Z  845  RAD 
Z  850  i 

Z855  i  routine  calculates  the  meAN  and  standard  deviation  of  an  array 
Z8G0  1  it  also  removes  the  mean  from  the  array 

Z8G5  i 

Z  870  Mean=0. 

Z  875  Sdv=0. 

Z880  FOR  1=1  TO  n 

Z885  Mean=Mean+Phird( I ) 

Z830  Sdv=Sdv+Phird( I )*Z 

Z89S  NEXT  I 

Z900  Mean=Mean/M 

Z905  Sdv=Sdv/M-MeanAZ 

Z910  IF  M<=  t  THEN  Sdv=( M/( M- 1 .  ) )*Sdv 

2915  FOR  1=1  TO  M 

Z9Z0  Phird( I )=Phird( I )-Mean 

Z9Z5  NEXT  I 

Z930  SUBEND 

Z935  i 

Z940  !  - 

Z945  ! 

Z950  SUB  Select <Limit,X,Y,Tf , Xrd< * ) , Yrd< * > , Trd< * ) , Urd< * > , Vrd( » > , Xopd< * > , Yopdt * ) 
, Topd( * ) ,Uopd( * ) ,Vopd( * ) , N , Nobs , Max_space_l ag , Max_t ime_l ag ) 

Z955  OPTION  BASE  1 
Z9B0  RAD 

Z965  ! 

Z970  !  -  routine  eliminates  the  distant  (in  space  and  time  >  points  and 

Z975  !  and  selects  the  most  "LIMIT"  near  points  to  an  interpolation 

Z980  !  point  X , Y , Tf . 

Z985  ! 

Z990  DIM  Index! 1089 ) ,Error( 1089) 

Z995  COM  /Ef ield/  El ,EZ 

3000  Cphase=0. 

3005  Nobs=0 

3010  R=FNUu(0.,0.  ,0. ) *FNVv<  0.  ,0. ,0. )-FNUv<0. ,0. ,0. >AZ 

3015  A1 1=FNVv<0. ,0. ,0. >/R 

30Z0  AZZ=FNUu(0. ,0. ,0. )/R 

30Z5  AIZ=FNUv(0. ,0. ,0. )/R 

3030  FOR  1=1  TO  N 

3035  Delt=Tf-Trd(  I ) 

3040  IF  A8S< Del t ) >Max_t ime_lag  THEN  3105 
3045  Delx=X-Xrd< I ) 

3050  Del y=Y-Yrd( I ) 

3055  R=SQR< ( Del x-Cphase*Del t )AZ+DelyAZ ) 

3060  IF  R>Max_space_l ag  THEN  3105 

3065  Nobs=Nobs+ 1 

3070  I ndex( Nobs  >  =  1 
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3075  C  -  I  =FNUuv  Del  a  ,  Del  y  .Del  L  ) 

3080  C.n 2  - FNUv(  Del  a  .  De  1  y  ,  De  1 1  ) 

3085  Cy!  '.  2 

3090  C  y  2  =  F  N  V  -v  (  De  1  a  .  Del  y,  Del  t  ) 

3095  Error'  Nobs>=< Cx 1' Z+Cy 1" Z ) * A  1 t  *Z . * ( Cx ! »CxZ+Cy 1 *CyZ  > *A 1 Z  +  <  CxZ “ 2 *CyZ ' Z ) *AZ Z 
3  1  00  Error  1  Nobs  i  =  ABS(  Z  .  -  Error!  Nobs  )  ) 

.3105  NEXT  I 

3110  IF  Nobs  - 0  THEN  3175 
3115  IF  Nobs rLim  t  THEN 

31Z0  CALL  Sort ! Error <*), I ndex (»>, Nobs > 

31Z5  Nobs=L  Lru.  t 

3130  END  IF 

31 35  FOR  1=1  TO  Nobs 

3140  J= Index! I) 

3145  Xopd( I >=Xrd( J) 

3150  Yopd! I )=Yrd( J ) 

3155  r0pd( I )  =  Trd<  J ) 

3160  Uopd<  I)=Urc*'  J) 

3165  Vopd(I)=Vrd<J> 

3170  NEXT  I 
3175  SUBEND 
3180  ! 

3185  !  - 

3190  ! 

3195  SUB  Sort! Cor! »), Index! *> ,N> 

3Z00  ! 

3Z05  1  -  A  shell  sort  routine  to  sort  index  and  cor  up  according  to 

3Z10  1  the  values  of  cor 

3Z15  ! 

3ZZ0  OPTION  BASE  1 
3ZZ5  RAD 
3Z30  Igap=N 

3Z35  IF  Igap<= 1  THEN  33Z0 

3Z40  Igap=Igap/Z 

3Z45  Imax=N-Igap 

3Z50  I ex=0 

3Z55  FOR  1=1  TO  Inax 

3Z60  Iplusg=I+Igap 

3Z65  IF  Cor! I )>=Cor( Iplusg)  THEN  3305 

3Z70  Save=Cor( I > 

3Z75  Cor! I )=Cor( Iplusg) 

3Z80  Cor! Iplusg )=Save 

3Z85  I save=I ndex! I ) 

3Z90  Index!  I )  =  Index< Iplusg) 

3Z95  Index! Iplusg)=Isave 

3300  I  ex  = 1 
3305  NEXT  I 

3310  IF  I  ex < >0  THEN  3Z50 
3315  GOTO  3Z35 
33Z0  SUBEND 
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33Z5  i 

3330  I - - 

3335  i 

3340  SUB  Veo t or_oa( Xopd( * ) , Yopd( * ) ,  Topd<  * ) ,Uopd( * )  ,Uopd( *  )  ,  N  ,  X  ,  Y , Tf ,  Uest .Vest  ,U 
err , Verr , I er ) 

3345  I 

3350  1  -  The  vector  space-tine  objective  analysis  routine. 

3355  i 

33G0  OPTION  BASE  1 
3365  RAD 

3370  DIM  A( 80,80) ,Cu( ZI78) ,  C v( 2 1 78  > 

3375  COM  /Ef ield/  El ,EZ 
3380  Uest=0. 

3385  Vest=0. 

3390  Uerr= 1 . 

3395  Verr=l. 

3400  I er=- 1 

3405  IF  NO0  THEN  GOTO  3615 

3410  CALL  Set_i nva(  A(  * ) , Xopd( * ) , Yopd( * > , Topd( * ) , N , I er ) 

3415  IF  Ier>0  THEN  3615 
34Z0  ! 

34Z5  !  -  calculate  the  correlation  matricies  CU  and  CV. 

3430  i 

3435  FOR  J=1  TO  N 
3440  Delt  =  Tf-Topd<  J ) 

3445  Delx=X-Xopd( J ) 

3450  Del y=Y- Yopd( J ) 

3455  Cu( J )=FNUu( Delx ,Dely ,Delt ) 

3460  Cu(  J+N  >  =  FNUv< Del  x  ,  Del  y  ,  Del  t ) 

3465  Cv< J+N)=FNVv(Deix,Dely,Delt ) 

3470  Cv( J)=Cu< J+N) 

3475  NEXT  J 
3480  i 

3485  I  -  calculate  the  error  variances  Uerr  and  Verr  and  velocity 

3490  !  estimates,  Uest  and  Vest. 

3495  i 
3500  NZ=N*Z 

3505  Uerr=0. 

3510  Verr=0. 

3515  FOR  1=1  TO  NZ 
35Z0  I pn=I +N 

35Z5  Eta=0. 

3530  FOR  J=1  TO  N 

3535  Jpn= J+N 

3540  Eta=Eta+A< I , J ) *Uopd( J ) +A( I , Jpn > * Vopdl J ) 

3545  IF  I <=N  THEN 

3550  Uerr=Uerr+Cu( I >  »A< I , J ) *Cu( J ) +Cu( Ipn)*A( I , Jpn)*Cv( J> 

3555  Verr=Verr+Cv< Ipn)*A( Ipn,Jpn)*Cv( Jpn)+Cu( Ipn)*A< I ,Jpn)*Cv( J) 

3560  ENO  IF 

3565  NEXT  J 
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3570  Uest=Uest  +  Cu< I  )*Eta 
3575  Vest  =  Vest  +  Cv( I  )»Eta 
'580  NEXT  I 

3585  Uest =  Uest +Uave 

3590  Vest=Vest+Vave 

3595  C a  au  =  FNUu( 0.  ,0.  ,0.  )+FNUv<0.  ,0. ,0.  ) 

3600  Cxx v=FNUv( 0. ,0. ,0. )+FNVv<0. ,0. ,0. ) 

3605  Uerr  = I . -Uerr/Cx au 
36)0  Verr-t . -Verr/Cxxv 
3615  SUBEND 
3620  ' 

3625  ! - 

3630  ! 

3635  SUB  Set_inva( A< * ) ,Xopd( * ) ,Yopd< * ) ,Topd< * ) .Nobs, Ier ) 

3640  ! 

3645  I  -  routine  sets  up  the  correlation  function  for  the  observations 

3650  i  given  the  positions  Xopd,  Yopd  and  times,  Topd,  it  returns  the 

3655  !  inverted  correlation  function  matrix. 

3660  ! 

3665  OPTION  BASE  1 

3670  RAD 

3675  DIM  Ip< 1089) 

3680  COM  /Ef ield/  El ,E2 

3685  Guard=1 .0E-30 
3690  D1 = 1 . 

3695  Ier=0 

3700  FOR  1=1  TO  Nobs 

3705  Ipnobs=I +Nobs 

3710  FOR  J=I  TO  Nobs 

3715  Jpnobs= J+Nobs 

3720  Delt  =  Topd( I >-Topd(  J) 

3725  Delx=Xopd(  I )-Xopd<  J ) 

3730  Dely=Yopd< I )-Yopd<  J) 

3735  f)<  I  ,  J  )=FNUu(  Del  x  ,  Del  y  ,  Del  t ) 

3740  A(  I  ,  Jpnobs  )=FNUv<  Delx  ,  Del  y  ,  Del  t ) 

3745  A( Ipnobs , J ) =A( I , Jpnobs ) 

3750  A( Ipnobs , Jpnobs ) =FNVv(  Del x , Dely , Del  t ) 

3755  NEXT  J 

3760  A<  I  ,  I  >=A( I , I >+E 1 

3765  A( Ipnobs .Ipnobs )=A< Ipnobs , Ipnobs  >+EZ 

3770  NEXT  I 
3775  NZ=Nobs*Z 
3780  i 

3785  1  -  invert  the  ( nobs*Z )*< nobs*Z )  matrix  A 

3790  i 

3795  CALL  Invmtx(  A<  *  )  ,  NZ  ,  A(  *  )  ,  NZ  ,  NZ  ,  D ,  I  p<  *  >  ,  I  er  ) 

3800  IF  D<Guard  THEN 

3805  PRINT  "Warning  the  determinant  is  very  small;  DET=  ";D 

3810  I er=- 1 

3815  END  IF 
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3820  SUBEND 
38Z5  i 

3830  i - 

3835  i 

3840  SUB  Est_mean( A( * ) ,Xopd( * ) ,  Yopd( * )  .  Topd( * ) ,N,U( *  > ,  V( * ) , Uave ,  Vave  > 
3845  ! 

3850  *  routine  calculates  the  estimated  mean,  and  then  removes 

3855  1  the  estimated  mean  from  the  observations. 

38G0  i 

38G5  OPTION  BOSE  1 
3870  RAD 

3875  DIM  C( Z 178) ,D<  Z 178) 

3880  CALL  Set_i nva( A< * ) , Xopd( * ) , Yopd( * ) , Topd( * > , N , I er ) 

3885  N2=N*2 

3890  FOR  1=1  TO  N2 
3895  I pn=I +P 

3900  C(I)=0. 

3905  D( I )=0. 

3910  FOR  K  =  1  TO  N 

3915  Kpn=K+N 

3920  C( I )=C( I )+A< I ,K)*U(K )+A( I ,Kpn>*V<  K ) 

3925  D( I >=D< I )+A(  I , K )+A(  I ,Kpn) 

3930  NEXT  K 

3935  NEXT  I 
3940  Sum1=0. 

3945  Sum2  =0. 

3950  Sum3=0. 

3955  Sum4=0. 

39G0  FOR  1=1  TO  N 
39G5  I pn=I +N 

3970  Suml =Sum1 +C< I ) 

3975  SumZ=SumZ+D< I ) 

3980  Sum3=Sum3+C( Ipn ) 

3985  Sum4=Sum4+D( I pn ) 

3990  NEXT  I 
3995  Uave=Sum1 /SumZ 
4000  Vave=Sum3/Sum4 
4005  ! 

4010  !  remove  the  calculated  means 

4015  I 

4020  FOR  1=1  TO  N 
40Z5  U( I )=U< I >-Uave 

4030  v< I >=V( I )-Vave 

4035  NEXT  I 
4040  SUBEND 
4045  ! 

4050  ' - 

4055  I 

40G0  SUB  Diag< Dpsi(  * ) ,M) 

4065  OPTION  BASE  1 
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407(5  RAD 
4075  i 

4080  IF  M  =0  THEN  4190 
4085  PRINT 

4090  PRINT  "Number  of  points:  "  ,M 
4095  Ave=0, 

4100  Sdv=0. 

4105  Psimax=Dp5i(  1  ) 

41  10  Psimin=Dpsi(  1  ) 

4115  Ave=Dps i <  1  ) 

4120  IF  M= 1  THEN  4155 
4125  FOR  1=2  TO  M 

4190  IF  Psimax<Dpsi(  I )  THEN  Psimax=Opsi(  I  ) 

4135  IF  Psimin>Opsi(  I )  THEN  Ps imin=Dps l < I ) 

4140  Sdv=<  <  I  -  2  )*Sdv+(  I  -  1  )*(  Dpsi<  I  )- Aver  >''2/1  )/(  I  -  1  )  . 

4145  Ave=( ( I- 1 )*AvetDpsi< I ) )/ I 

4150  NEXT  I 

4155  Rms=SQR( Sdv+Ave' 2 ) 

4160  PRINT  "Mean  and  Variance:  ",Ave,Sdv 
4165  Sdv=SQR( Sdv ) 

4170  PRINT  "Standard  Deviation:  “,Sdv 

4175  PRINT  "RMS  of  field:  ",Rns 

4180  PRINT  "minimum  of  field:  " .Psimin 

4185  PRINT  "maximum  of  field:  “.Psimax 

4190  SUBEND 
4195  1 

4200  !  - - 

4205  ! 

4210  SUB  Invmtx( A( * ) ,Na,V(  * ) , Nv ,N,D, Ip(  * ) , ler ) 

42  15  l 

4220  !  -  routine  inverts  the  Matrix  A. 

4225  ! 

4230  OPTION  BASE  1 

4235  RAD 

4240  Iexmax=75 

4245  Ier=FNIerinv(  N,Na,Nv) 

4250  IF  lcr<>0  THEN  4600 
4255  FOR  J=l  TO  N 
4260  Ip(J)=0 

4265  FOR  1=1  TO  N 

4270  V< I  ,  J)=A( I , J) 

4275  NEXT  I 

4280  NEXT  J 
4285  D= 1 . 

4290  lex=0 

4295  FOR  M= 1  TO  N 
4300  Vmax=0. 

4305  FOR  J=l  TO  N 

4310  IF  Ip(J)<>0  THEN  4355 

4315  FOR  1=1  TO  N 
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4320  IF  Ip( I  )  v  0  THEN  4350 

4325  Vh=ABS< V< I , J> > 

4330  IF  Vnax'Vh  THEN  4350 

4335  Vrcax=Vb 

4340  K=I 

4345  L=J 

4350  NEXT  I 

4355  NEXT  J 

4350  I  p(  L  )  =K 

4355  Npm=N+M 

4370  I p( Npm ) =L 

4375  D=D*y <  K , L ) 

4380  IF  ABS<  D) <= 1 . 0  THEN  4400 

4385  D=D* . 1 

4390  Iex=Iex+1 

4395  GOTO  4380 

4400  Pvt=V( K , L ) 

4405  IF  M= I  THEN  Pvtmx=ABS< Pvt ) 

1410  IF  < ABS< Pvt/M)+Pvtmx >=Pvt*x  THEN  4570 

4415  V(  K  , L  >  =  1  . 

4420  FOR  J=1  TO  N 

4425  Hoi d=V< K , J ) 

4430  V(K,J)=V(L,J) 

4435  L , J)=Hold/Pvt 

4440  NEXT  J 

4445  FOR  1=1  TO  N 

4450  IF  I =L  THEN  4480 

4455  Hoi d=V( I , L ) 

4450  V(I,l)=0. 

4455  FOR  J=l  TO  N 

4470  V< I  ,  J ) =V< I , J)-V(L,J)*Hold 

4475  NEXT  J 

4480  NEXT  I 

4485  NEXT  M 

4490  M=N+N+ 1 

4495  FOR  J=1  TO  N 

4500  M=M- 1 

4505  L=Ip( M) 

4510  K=Ip( L ) 

4515  IF  K=L  THEN  4550 

4520  D=-D 

4525  FOR  1=1  TO  N 

4530  Hoi d=V( I , L ) 

4535  V( I ,L  )=V< I ,K  > 

4540  V( I ,K )=Hold 

4545  NEXT  I 

4550  NEXT  J 

4555  IF  Iex>Iexniax  THEN  4585 
45G0  0=0* 10. A I  ex 

4565  GOTO  4600 
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457®  Ier=33 

4575  PRINT  "MATRIX  SINGULAR  IN  INVMTX;  ERROR  CODE  IS  " .Ier 
4580  GOTO  4S00 
4585  I er  ~ I 
4590  D= I  ex 

4595  PRINT  "DETERMINANT  TOO  LARGE  IN  INVMTX;  ERROR  CODE  IS  ",Ier 
4G00  SUBEND 
4G05  i 

4G10  i  - 

4615  i 

46Z0  DEE  FNIennv(  N.Na.Nv) 

46Z5  OPTION  BASE  I 

4630  RAD 

4635  Iennv=0 

4640  IF  N>= 1  THEN  4660 

4645  Ierinv=34 

4650  PRINT  "N  <  1  IN  INVMTX" 

4655  GOTO  4695 

4660  IF  Na>  =  N  THEN  4680 

4665  Iennv=35 

4670  PRINT  " NA<N  IN  INVMTX" 

4675  GOTO  4695 

4680  IF  Nv>=N  THEN  4695 

4685  Iennv=36 

4690  PRINT  "NV  <  N  IN  INVMTX" 

4695  RETURN  Iermv 

4700  FNEND 
4705  ' 

4710  ! 

4715  i 

47Z0  Get_data: SUB  Get_data<  F  x 1 e name$ , Xdata< »  > , Ydatal  *  > , Tdata( * ) , Udata( *  > , Vdata( * 
) , Ndaia ) 

47Z5  ! 

4730  1  -  routine  reads  in  the  observations. 

4735  ' 

4740  ASSIGN  @Path_in  TO  Filenames 
4745  Ndata=0 

4750  PRINT  "  **********  Echo  check  of  the  first  ten  records  »***###*** 

*  " 

4755  PRINT  "  No.  X_POS  Y_P0S  TIME  U-COMP. 

V-COMP.” 

4760  ON  END  @Path_in  GOTO  4815 

4765  DISP  USING  "K ,K , K" ; "Reading  the  observed  data  from  ".Filenames,",  please 
wait . " 

4766  LOOP 

4770  ENTER  @Path_in; X , Y , T , U , V 

4775  IF  Ndata<= 1 0  THEN  PRINT  USING  "DODD, X , MD. DDDDDDE , X , MD. DDDODDE , X , MD. DODDDD 
□DDE , X , MD. DDDDDDE , X , MD. DDDDDDE " ; Ndat a , X , Y , T , U , V 
4780  Ndata=Ndata+ 1 
4785  Xdata<  Ndata)=X 
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4790  Ydata( Ndata) =Y 
4795  T  dat a( Ndat  a )  =  T 
4800  Udata( Ndata)=U 
4805  Vdata( Ndata ) = V 
4810  END  LOOP 

4815  ASSIGN  @Path_in  TO  *1  Closing  input_file. 

4820  DISP 
4825  SUBEND 
4830  ! 

4835  i 
4840  I 

4845  Get_data:SUB  Get_ip( F i 1 ename® , Xdat a< * > , Ydat a( * ) , Ni p ) 

4850  i 

4855  !  - routine  reads  in  the  int er/ex t rapol at  ion  positions. 

48G0  I 

48G5  ASSIGN  @Path_in  TO  Filename® 

4870  PRINT  "*  Echo  check  of  the  first  ten  records  »" 

4875  PRINT  "  No.  X_IP_POS.  Y_IP_P0S." 

4880  Nzp=0 

4885  ON  END  @Path_in  GOTO  Close_file 

4890  DISP  USING  " K , K , K " j "Reading  the  interpolated  positions  from  ".Filename®," 
,  please  wait." 

4891  LOOP 

4895  ENTER  @Path_in;X,Y 

4900  IF  Nip<= 1 0  THEN  PRINT  USING  " DDOD . XX , MD. DDDDDDE , XX . MD. DDDDDOE"  ;  Nip , X , Y 

4905  Nip=Nip+1 

4910  Xdata(Nip)=X 

4915  Ydata( Nip ) = Y 

4920  END  LOOP 

4925  Close_f ile: ASSIGN  @Path_in  TO  *  !  Closing  input_file. 

4930  DISP 
4935  SUBEND 
4940  i 
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10  i  PROGRAM  SOAOI FF 

!  1  ' 

20  1  Program  Calculates  Oiference  Between  Two  15x15  Arrays 

30  ' 

1000  OPTION  BASE  1 

1 00B  DIM  F i 1 e$l 40] , Fi t 1 e 1 $1 801 ,TitleZ$[80] 

1010  REAL  Farea 1 (  1 5 , 1 5 ) , FareaZ (  1 5 . I  5 ) , Farea(  I  5 , I  5 ) , F area_sum 
1015  i 

1020  OUTPUT  KBD; CHR$<  Z55 )&CHR$( 75  ) ; 

1025  F l 1 e$= " “ 

1030  LINPUT  "Enter  first  file  to  compute  difference  (include  device)  exit  ",F 

1 1  e® 

1035  F i 1 e$=  TR! M$( UPC$<  F i 1 e$  > ) 

1040  IF  File$=""  THEN  GOTO  Stopprogram 
1045  ON  ERROR  GOTO  Data_error 
1050  ASSIGN  ©File  TO  File® 

1055  ENTER  @F i 1 e ; T 1 1 1  e 1  $ , T i 1 1  e2$ , Farea I  (  * ) 

1060  ASSIGN  ©File  TO  * 

1065  ' 

1070  F l 1 e$= 

1075  LINPUT  "Enter  second  file  to  compute  difference  (include  device)  <exit>". 
File® 

1080  F ile®=TRIMS< UPC$< File®)  > 

1085  IF  File$=""  THEN  GOTO  Stop_program 
1080  ASSIGN  ©File  TO  File® 

1095  ENTER  @File;TitlelS,TitleZ$,Farea2(  * ) 

1 100  ASSIGN  ©File  TO  * 

1110  F i I e$= " ” 

1115  LINPUT  "Enter  file  for  File!  -  FileZ  difference  (include  device)  <exit>", 
File® 

11Z0  F i le*=TRIM$( UPCS( File®  )  ) 

1tZ5  IF  File®=""  THEN  GOTO  Stop_program 

1130  ON  ERROR  GOTO  1 140 

1135  CREATE  BDAT  File®, 10 

1140  ON  ERROR  GOTO  Data_error 

1145  ASSIGN  ©File  TO  File® 

1150  MAT  Farea=  Fareal -FareaZ 

1155  OUTPUT  ©FilesTitlelS.TitleZS, Fareal  * ) 

1 160  ASSIGN  ©File  TO  * 

1165  MAT  Farea=  Farea  .  Farea 
1170  Farea_sum=SQR(  SUM< Farea) ) 

1175  DISP  "Curomulated  Dif f erence : " ; Farea_sum 
I  180  ASSIGN  ©File  TO  * 

1185  STOP 

1190  Stop_program:  1 
1195  DISP  "Stopped" 

1Z00  STOP 

IZ05  Oata_error:  ! 

1Z10  DISP  ERRMS) " ;  File  ’ " ; F i le$; " ’ " 

1 Z 1 5  END 
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1000  i  PROGRAM  PLOTJJA 
1005  i 

1010  i  -  Procran  plots  the  objective  analysis  field. 

1015  ' 

1 0Z0  OPTION  BASE  1 

1 02 5  DIM  T i 1 1 e I $[ 80]  ,TitleZ$I80]  ,User$I40]  ,Ansu$( 10]  ,Old$140] , Oa_f 1 1 e$t 40] 

1030  DIM  Sfc < 15,15) 

1035  i  - 

1040  i  Assign  the  Objective  analysis  results  file  to  Oa_file$ 

1045  0 1 d$= " " 

1050  LOOP 

1055  OUTPUT  KBD; CHR$(  2 55 ) &CHR$( 75 ) ; 

1 060  LINPUT  "Enter  the  O.A.  data  file  to  be  plotted  <ex 1 1 > " , Oa_f i 1 e$ 

1 0G5  Oa_f i 1 e$  =  TRI M$(  UPC$<  Oa_f ile®> ) 

1070  EXIT  IF  Oa_f 1 1  e$=" " 

1075  i  - 

1080  1  Input  the  x-  and  y-gnd  points  of  conputating  field. 

1085  INPUT  "Enter  the  x-  and  y-grid  points  of  conputating  field,  (Imax.Jmax)" 
I  max , Jmax 

1090  IF  I  max  <  -  Z  1  OR  JmaxOZI  THEN  It  30 
1095  BEEP 

1100  PRINT  TABXY( 10,10),"  ***  Warning  ***" 

1105  PRINT  TABXY( 1 0 , 1 1 ) , " The  sizes  of  array  Sfc(i,j)  defined  in  this  version" 

1110  PRINT  TABXY( 1 0, 1 Z ) , " is  incorrect,  make  correction  in  lines  1015  and  1055 

1115  PRINT  TABXY(  1 0 , 1 3 ), "Program  execution  is  halted." 

1 1 20  STOP 

1125  1  Set  the  minimum  and  maximum  contours 

1130  INPUT  "Enter  the  minimum  contour  to  be  plotted. " ,Min 

1135  INPUT  "Enter  the  Maximum  contour  to  be  plotted. ” .Max 

1140  i  - - - 

1145  l  Set  number  of  contours(Nc)  or  set  contour  interval ( C l ) ,  if  Nc=0 

1150  INPUT  "Enter  the  number  of  contour  or  the  Contour  interval", Nc 

1155  IF  INT(  Nc* 100/ 100>=Nc  THEN 

1160  C i =ABS( Max -Mi n ) /Nc 

1165  ELSE 

1170  C i =Nc 

1  175  END  IF 

1180  !  - 

1185  1  user  1 abel 

1190  User S=" " 

1195  LINPUT  "Enter  the  user's  label  <40  char  max)", User® 

1 Z00  UserS=TRIM$< User®) 

1 Z05  i  - 

1210  1  screen  or  plotter 

1215  Answ$="S" 

1ZZ0  INPUT  "Plot  on  screen  or  pen  plotter  (S/P)  (screen)" , Ansu® 

1 ZZ5  Dvc=0 

1 Z 30  IF  UPCS< Ansu®[ 1 ; 1 1 >  =  "P"  THEN  Dvc=1 

1 Z  35  I  - 
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1Z40  1  dump  graph  ics  to  printer 

1Z45  Oump_ graph-0 
1 Z50  IF  Dv.;=0  THEN 
I  Z 55  Ar.su*-- "N" 

IZS0  INPUi  "Dump  the  graphics  to  printer  (Y/N)  no>",Ansu$ 

1 Z 65  IF  UPC$< AnsuSI  1  ;  I  1  )  =  " Y"  THEN  Oump_graph=l 

1Z70  END  IF 

I Z75  i  - 

1 Z  80  C$=CHR$< Z5& ) &" K " 

I  Z  85  i  - 

1Z90  1  read  in  data 

1  Z 95  IF  01d$00a_f  lie*  THEN 

1300  OUTPUT  Z  USING  "#,K"iC$  i  Clean  the  screen 

1305  ASSIGN  @Path_in  TO  Oa_file* 

1310  ON  END  @Path_in  GOTO  1505 

1315  ENTER  @Path_i n ; T i 1 1 e 1  $ 

1  3Z0  ENTER  @Path_in; TitleZ* 

1 3Z5  FOR  I  =  I  max  TO  I  STEP  -1 
1330  FOR  J=l  TO  Jmax 

1335  ON  END  @Path_in  GOTO  1350 

1340  ENTER  @Path_in;Sfc< I , J) 

1345  NEXT  J 

1350  NEXT  I 

1355  OFF  END  @Path_in 

1350  ASSIGN  @Path_in  TO  * 

1355  0'ld$=0a_f  ile* 

1370  END  IF 

1 375  OUTPUT  Z  USING  "#,K";C$  !  Clean  the  screen 

1380  IF  Dvc=0  THEN 

1385  GINIT 

1330  GCLEAR 

1395  GRAPHICS  ON 

1400  ELSE 

1405  PLOTTER  IS  705,"HPGL" 

1410  END  IF 

1415  Contour! Sfc( *) ,Min .Max ,Ci , 1 ., Title  1 $, Tit leZ-  U  ')!  do  the  plotting 

1 4Z0  BEEP 

14Z5  IF  Ovc-0  THEN 

1430  ON  KBD  ALL  GOTO  1440 

1435  Idle:GOTO  Idle  I  vieu  the  plot  as  long  as  you  want. 

1440  OFF  KBD 

1445  IF  Dump_graph  THEN 

1450  GRAPHICS  OFF 

1455  DISP  "Dumping  graphics  ..." 

1450  DUMP  GRAPHICS  CRT  TO  ttPRT 

1465  END  IF 

1470  GCLEAR 

1475  GRAPHICS  OFF 

1480  END  IF 

1485  END  LOOP 
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1490  DISP  "Stopped” 

1495  STOP 

1 500  i  - 

1505  OFF  END  @Path_in 
1510  ASSIGN  @Path_i n  TO  * 

1515  DISP  ERRM$; ” ;  File  ’ " ; Oa_f ile$; " ’ " 

I 5Z0  END 
1 525  i 
15.30  i 
1535  i 

1 540  Contour: SUB  Contour ( Sfc( * ) ,  Mi n , Max .Interval .Extreme s,Titlel$,TitleZ$, Us er$> 
1545  OPTION  BASE  1 
1550 
1555 
1  5G0 
15G5 
1570 
1575 
1580 
1585 
1590 
1595 
1B00 
1605 
1610 
1615 
1620 
1625 
1630 
1635 
1640 
1645 
1650 
1655 
1660 

1665  INTEGER  I , J , I  max , Jmax 

1670  COM  /G_units/  Gdu_xmax , Gdu_ymax , Udu_xmi n , Udu_xmax , Udu_ymin , Udu_ymax , Show 
1675  CALL  6du( X_gdu_max , Y_gdu_max , Xmid.Ymid) 

1680  VIEWPORT  0. , X_gdu_max ,0, Y_gdu_max 
1685  WINDOW  0. ,X_gdu_max ,0,Y_gdu_max 
1690  LINE  TYPE  1 

1695  CALL  Label ( 4 , .6 ,0 ,5 , 1  , Xmid , . 98*Y_gdu_max , Tit le 1  $ ) 

1700  CALL  Label ( 4 , . 6 ,0 . 5 , I  , Xmid , . 94* Y_gdu_max , Tit  1 eZ$ ) 

1705  Min_$=VAL$(  INT(  MINI Sf c( * ) )* 1 000 > / 1 000 ) 

1710  Max_$=VAL$< INT< MAX< Sf c<  * ) > *  1 000 ) / 1 000 > 

1715  C i_$=VAL$(  INTI  I nterval *  1 000 ) / 1 000 > 

17Z0  CALL  Label < 3 . 5 , . 6 ,0 , 5 , 1 , Xmid , . 1  * Y_gdu_max , "Min. :  "&Min_$&";  Max.:  “&Max 
$&";  Cl:  "&Ci_$> 

17Z5  CALL  Label! 3.5. .6.0.5. 1 .Xmid. .05»Y  qdu  max.User$) 

1730  VIEWPORT  0. , X_gdu_max , . 1 4* Y_gdu_max , Y_gdu_max* . 9 


Copyright  1983.  Hewlett-Packard  Company 
All  Rights  Reserved 

This  subprogram  plots  a  contour  map  of  the  array  Sfc(»),  and 

optionally  plots  local  minima,  maxima,  and  statistics. 

Sfc(*):  This  is  the  two-dimensional  real  array  containing  the 

data  to  be  plotted.  It  need  not  be  square. 

Min  (i  Max:  These  are  the  lowest  and  highest  levels,  respect  i  vel  y , 

of  the  contour  lines.  These  allow  you  to  specify  the 
exact  range  within  which  you  want  contours.  Every 
contour  line  outside  of  this  range  will  not  be  plotted. 

Interval:  This  specifies  how  far  apart  the  contour  lines  have  to 

be  (in  value,  not  in  distance).  The  smaller  the  inter¬ 
val,  the  denser  the  contour  plot. 

Extremes:  This  is  a  logical  variable  which  specifies  whether  or 

not  to  label  local  maxima  and  minima.  A  local  maximum 
is  a  point  whose  value  is  larger  that  its  eight 
neighbors  immediately  to  the  west,  northwest,  north, 
northeast,  east,  southeast,  south,  and  southwest.  A 
local  minimum  has  a  corresponding  definition. 
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1735  Imax=5IZE( Sfc  ,  I  ) 

1740  Jma*=SI  ZE(  Sf c  .  Z  ) 

1745  Shouf  1  , ( Jmax  > . <  I  max ) , 1  ) 

1750  LINE  TYPE  i 

1755  PEN  1 

1760  CLIP  1 . Jmax , Imax , 1 
1765  AXES  I  ,  I  ,  I  ,  1  ,  I  .  1  .  I 
1770  AXES  1  . I  .Jmax ,  Imax , 1 . 1  ,  1 
1775  Northeast=0  1  \ 

1780  Northuest=1  !  >  Figure  what  to  do  for  case  4. 

1785  Cross=0  '  / 

1790  PEN  Z 

1795  FOR  1=1  TO  Imax-1 

1800  FOR  J=1  TO  Jmax-1 

1805  Big=MAX( Sf  c( I , J ) , Sf c< I , J+ 1 )  ,Sfc( 1  +  1 , J) ,Sfc< 1+1 ,J+1 > ) 

1810  Small =MIN< Sfc< I , J) ,Sfc< I , J+ 1 > , Sf c< I  +  1  , J) ,Sfc< 1  +  1 , J+1  )  ) 

1815  FOR  Cont=Min  TO  Max  STEP  Interval 

1 8Z0  IF  Cont  >Smal 1  AND  Cont<Big  THEN 

1 8Z5  LINE  TYPE  1 

1830  IF  Cont <0  THEN  LINE  TYPE  4 

1835  Top=Cont>MIN( Sfc< I , J) ,Sfc< I , J+1 ) )  AND  Cont<MAX< Sfc< I , J) ,Sfc( I , J+1 > ) 

1840  Bottom=Cont>MIN( Sfc< 1  +  1  ,J)  ,Sfc< 1  +  1 , J+1  )  )  AND  Cont <MAX< Sf c< I + 1 , J ) , Sf c< 

1+1 , J+1 ) ) 

1845  Lef t  =  Cont >MIN(Sfc( I ,J) ,Sfc< 1+1 ,J) )  AND  Cont <MAX( Sf c( I , J ) , Sf c< I  + 1 , J > ) 

1850  Right  =  Cont>MIN< Sfc< I , J+1 ) ,Sfc< 1  +  1  , J+1  )  )  AND  Cont <MAX( Sf c( I . J+ 1 ) , Sf c< I 

+1 , J+1 ) ) 

1855  SELECT  Top+Bottom+Lef t+Right 

1860  CASE  0  !  Do  nothing . 

1865  CASE  Z  1  Tuo  intersections,  so  draw  one  line . 

1870  IF  Top  THEN 

1875  Jtop=J+( Cont-Sf c(  I ,J> )/(Sfc< I , J+1  >-Sfc< I , J) ) 

1880  IF  Bottom  THEN  !  Top  and  Bottom . 

1885  Jbottom=J  +  (Cont-Sfc< I  + 1 , J > )/<  Sf c< 1  + 1 , J+ 1 >-Sfc< 1  +  1  ,  J)  ) 

1890  MOVE  Jtop.I 

1895  DRAW  Jbottom ,1+1 

1900  ELSE  !  (not  Bottom) 

1905  IF  Left  THEN  !  Top  and  Left . 

1910  Ileft=I+(Cont-Sfc< I , J) )/( Sfc( 1  +  1 ,J)-Sfc( I , J) ) 

1915  MOVE  Jtop.I 

1920  DRAW  J , 1 1 ef t 

19Z5  ELSE  !  Not  left,  therefore  Top  and  Right . 

1930  Ir ight  =  I  +  ( Cont-Sf c( I ,J+1 >  > / (  Sfc(  1+1 ,J+1 )-Sfc(  I ,J+1  ) ) 

1935  MOVE  Jtop.I 

1940  DRAW  J+1 , Iright 

1945  END  IF  !  (  if  left ) 

1950  END  IF  !  ( if  bottom) 

1955  ELSE  !  (not  Top) 

1960  IF  Bottom  THEN 

1 965  Jbot tom= J  +  (  Cont-Sfc< 1  +  1 ,J> )/(Sfc(  1  +  1  ,J+1  )-Sfc(  1  +  1  ,J) ) 

1970  IF  Left  THEN  !  Bottom  and  Left . 
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1975 
1980 
1985 
1  990 

1  995 
2000 
2005 
2010 
2015 
2020 
2025 
2030 
2035 
2040 
2045 
2050 
2055 
2060 
2065 
2070 
2075 
2080 
2085 
2090 
2095 
2100 

2  105 
2110 
2115 
2120 
2125 
2  130 
2135 
2140 
2145 
2150 
2155 
2160 
2  165 
2170 
2175 
2  180 
2185 
2  190 
2195 
2200 
2205 
2210 
2215 
22Z0 


I  lef  t=I  +(  Cont-Sf  c<  I  ,  J  >  >  /  (  Sf  c(  I  +  1  ,  J  )-Sf  <:(  I  ,  J  )  > 

MOVE  J, lief t 
DRAU  Jbottom,  1  +  1 

ELSE  1  Not  left,  therefore  Bottom  rand  Right. 

I r lght  =  I  +  ( Cont-Sf c< I ,J+1  ))/<Sfc( I  +  l  ,J+I  )-Sf c( I , J+ 1  >> 
MOVE  Jbottom ,1+1 
DRAW  J+ I . I r lght 
END  IF  i  <  if  lef t ) 

ELSE  1  Not  Bottom,  therefore  Left  and  Right... 

I lef t  =  I  +  ( Cont-Sf c( I  ,  J  )  )/( Sf c( I  +  1  , J )-Sf c( I , J ) ) 

I r ight  =  I +(  Cont-Sfc< I , J+ 1 ) )/( Sf c< I  +  1 , J+ 1 ) - Sf c< I  ,  J+ 1  > ) 
MOVE  J.Ileft 
DRAW  J+ I , Ir ight 
END  IF  I  ( if  bottom) 

END  IF  !  ( if  top) 

CASE  4  1  Four  intersections . 

Jtop=J+<  Cont-Sf c(  I , J) )/<  Sfc( I , J+1 )-Sfc( I , J) ) 

Jbot  tom= J+(  Cont-Sf c( 1  +  1 , J ) )/<  Sf c< 1  + 1 ,J+1 )-Sfc<I  +  1 ,J>> 

I lef  t  =  I  +  ( Cont-Sf  c( I  ,  J ) >/< Sfc< 1  +  1 , J)-Sfc( I , J) ) 

Iright=I+<  Cont-Sfc< I , J+1 ) )/< Sfc( 1  +  1  , J+1 )-Sfc( I  ,  J+1  )  > 

IF  Northeast  THEN 
MOVE  J.Ileft 
DRAU  Jt  op , I 
MOVE  Jbottom,I+1 
DRAU  J+1 , 1  right 

END  IF  !  < if  northeast) 

IF  Northuest  THEN 
MOVE  J.Ileft 
DRAU  Jbottom ,1+1 
MOVE  Jtop.I 
DRAU  J+1  , 1  right 

END  IF  1  < if  northwest) 

IF  Cross  THEN 
MOVE  J.Ileft 
DRAU  J+1 , 1  right 
MOVE  Jtop.I 
DRAU  Jbottom, I+l 
END  IF  I  < if  cross) 

END  SELECT 
END  IF 
NEXT  Cont 
NEXT  J 

NEXT  I 

IF  Extremes>0  THEN 

I 

LINE  TYPE  1 
Image$="K" 

FOR  1  =  2  TO  I  max- 1 
FOR  J=2  TO  Jmax-1 

Point=Sf c( I , J )  i  (The  point  we’re  working  on) 


B-66 


n  ->  "t  c 


Min=MIN(  Sf<:(  I  -  I  ,  J  I  )  ,  Sf  c(  I  -  I  ,  J  >  ,Sfc<  1-1  .  J  +  1  )  ,Sfc<  I  ,  .J-  1  )  ,Sf  c<  I  ,JH  )  .  ;i  F  . . 
< 1  +  1  , J  -  1  >  .Sfc( I  +  1  . J )  Sfc< I H  . J+ I  ) ) 

22  30  Max  =MA,X( Sfc(  f- 1  , J~  I  ) ,Sfc<  [ -  I  ,  J)  ,Sfc( I  -  I  ,  J+ I  )  , Sfc(  I  ,  J- 1  )  ,  Sfc<  I , J+  I  )  , Sf<; 

(1  +  1  .  J-  1  )  ,Sf c(  1  +  1  .  J  )  . S f  c (  1  +  1  , J+ 1  )  ) 

Z  Z  35  IF  Point  Max  OR  Point. Min  THEN 

2240  CALL  Label( 1 ,.6,0,5,4,(J),( I ),"+") 

2245  [F  Point: Max  THEN 

2ZS0  CALL  LabeK  3,  .6, 0.5. 4, <  J)  ,(  I  )  ,"H"  ) 

2  255  ELSE 

2260  CALL  Label! 3, .G,0,5,4.< J) .< I ) ,"L" ) 

2  265  END  IF 

ZZ70  IF  Extrenes>l  THEN 

2711  CALL  LabeK  3  , .  G  ,  0 , 5 , 4  ,  J  ,  I  +  .  Z  , )  1  No  label,  just  setup 

2280  LABEL  USING  I nage$ ; Poi nt 

2285  END  IF 

2290  END  IF 

2295  NEXT  J 

2300  NEXT  I 

2305  END  IF  i  (if  Stats ) 

Z310  MOVE  X_gdu_max , V_gdu_max 
2315  PEN  0 
2320  SUBEND 
2325  i 

2330  |  - 

2335  i 

2340  Shou:SUB  Shou< Xlef t , Xr ight , Y1 ou , Yhigh > 

Z345  I  This  simulates  the  system  command  SHOW,  but  saves  the  variables  so 

2350  1  the  routines  Setgu  and  setuu  work 

2355  COM  /G_umts/  Gdu_xmax ,Gdu_ymax , Udu_xmin , Udu_xmax , Udu_ymin , Udu_ymax , Shou 

Z 360  IF  Gdu_xmax=0  THEN 

Z3G5  Gdu_xmax=100*MAX<  1 .RATIO) 

2  370  Gdu_ymax= 1 00*MAX(  1  , 1 /RATIO) 

2375  END  IF 
2380  Udu_xmin=Xlef t 

2385  Udu_xmax=Xr ight 
Z390  Udu_ymin=Ylou 
Z395  Udu_ymax=Yhigh 
2400  Shouj=  I 

2405  SHOU  Xleft  .Xright  .YIouj,  Yhigh 

Z410  SUBEND 

2415  Gdu:SU8  Gdu( X_gdu_max . Y_gdu_max , OPTI ONAL  Gdu_xmid , Gdu_ymid ) 

2420  1  This  returns  Xright,  Yhigh  and  their  respective  midpoints  in  GDUs. 

Z4Z5  1  Note  that  if  Gud_xmid  is  defined,  Gdu_ymid  must  be  also. 

2430  COM  /G_units/  Gdu_xmax , Gdu_ymax , Udu_xmin , Udu_xmax , Udu_ymi n , Udu_ymax , Show 

2435  IF  Gdu_xmax=0  THEN 

2440  Gdu_xmax=100*MAX( 1 .RATIO) 

2445  Gdu_ymax  =  100«MAX< 1 , 1/RATI  0) 

2450  END  IF 

2455  X_gdu_max=Gdu_xmax 

Z4B0  Y_gdu_max=Gdu_ymax 
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7465  [F  NPAR2  THEN 

7470  Gdu_xmid=Gdu_xmax* . 5 

7475  Gdu_ym id=Gdu_ymax * . 5 

7480  END  IF 

2485  SUBEND 

Z490  i 

2495  ! 

7  500  ' 

Z505  Label:SUB  Label ( Cs lze , Asp_rat 10 , Ldir , Lorg , Pen , X , Y , Text® ) 

2510  1  This  defines  several  systems  variables  (in  CSIZE,  LDIR,  etc.),  and 
7515  1  labels  the  test  (if  any)  accordingly. 

7520  DEG 

2525  CSIZE  Cs l ze , Asp_r at  10 

2530  LDIR  Ldir 

2535  LORG  Lorg 

2540  PEN  Pen 

2545  IF  Text$<>"“  THEN 

2550  MOVE  X , Y 

2555  LABEL  USING  "#,K"iText$ 

2560  END  IF 

2565  PENUP 

2570  SUBEND 
2575  i 
2580  ! 

2585  ! 

2590  Scale:SUB  Scale( Surface! *> , Neu_min,New_max > 

2595  OPTION  BASE  1 

2600  !  This  routine  scales  a  matrix  such  that  it  will  have  a  new  lowest 
2605  1  value  of  Neu_min  and  a  new  higest  value  of  New_max . 

2610  DISP  USING  "K'VScaling  the  surface  array  from  "  ,  Neu_min ,  "  to  "  ,  Neu_max  ,  " 

2615  Min=MIN( Surf ace( * ) ) 

2620  Max=MAX< Surf ace( * ) ) 

2625  IF  Min=Max  THEN  1  Array  is  completely  flat 
2630  MAT  Surface=  <Neu_min) 

2635  SUBEXIT 

2640  END  IF 

2645  MAT  Surface=  Surface-! Min) 

2650  Range_recip=(  New_max-Neu_min)/(  Max-Min) 

2655  MAT  Surface=  Surf ace*< Range_recip > 

2660  MAT  Surface=  Surf ace+( New_mi n ) 

2665  DISP 

2670  SUBENO 
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